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General Notation 


n-dimenslonal state vector of discrete linear stochastic 
process 


m-dimensional observation vector 


k-difflensional zero mean white Gaussian process noise 


in-dlmensional zero mean white Gaussian data noise, Indepen* 
dent of w(i) 

nxn state transition matrix 


aT(i) 


nxk noise transition matrix 


mxn data coefficient matrix 


n-dimensional vector of data coefficients 


k-dimensional colored process noise vector 
kxk diagonal colored noise transition matrix 
b-dijuensional vector of bias parameters 
mxm diagonal covariance of data errors, v(i) 


covariance of scalar data error 


kxk diagonal covariance of white process noise 


minimum variance estimate of x(i) given data {z(J)}, J 1 i 

a priori estimate of x(i), based upon data {z(J)), j < i 

nxm filter gain matrix 

nxn error covariance matrix for 8(i) 

nxn error covariance matrix for St(i) 


innovations error covariance 


nxn square root of estimate error covariance matrix 
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Statistical and Filter-Related Quantities (contd) 


U nxn unit upper triangular factor of error covariance matrix 

D nxn diagonal matrix, related to error covariance matrix 

by the identity P s ODU” where U is defined above 


Ujj, Djj those portions of the U-D covariance factors which correspond 
to the X parameters 


0 ^ standard deviation ot tno parameter x 

E{*} expectation operator 

X - N(0,F) X is a normally distributed random variable with zero mean 
and covariance P 


Mathematical Symbols 
e is an element of 

the set consisting of n- tuples whose components are real 
numbers 

1 summation 

r 

= rounded to 

^ not equal to 

1 1 ‘ 1 1 Euclidian vector norm 

1 1*1 Ip with metric D, i.e., I|v|lp s^v^Dv^^^^ 

( sequence of n quantities 

:s replace in computer storage 

>> much greater than 

Abbreviations 

WGS weighted Gram-Schmldt 

MWGS modified weighted Graro-Schmidt 

RSS root sum of squares 

GS Gram>Schmldt 
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km/ sec 


kilometer per second 


m/sec 
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gravitational mass 


Superscripts 


matrix transpose 


matrix inverse 


upper triangular Cholesky square root 

inverse of upper triangular Cholesky square root 

after incorporating measurement 

before incorporating measurement 

after jth iteration (section 3*3) 

after (n-j)th iteration (sections 3.^ and 3-5) 
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Abstract 

In this report an improved computational form of the discrete 
Kalman filter is derived using an upper triangular factorization of 
the error covariance matrix. The covariance P is factored such that 
P = UDU^ frtiere U is unit upper triangular and D is diagonal. Recursions 
are developed for propagating the (J-D covariance factors together with 
the corresponding state estimate. The resulting algorithm, referred 
to as the U-D filter, combines the superior numerical precision of 
square root filtering techniques with an efficiency comparable to that 
of Kalman's original formula. Moreover, this method is easily implemented 
and involves no more computer storage than the Kalman algorithm. These 
characteristics make the U-D method an attractive real-time filtering 
technique . 

A new covariance error analysis technique is obtained from an 
extension of the U-D filter equations. This evaluation method is flex- 
ible and efficient and may provide significantly improved numerical 
results. Cost comparisons show that for a large class of problems 
the U-D evaluation algorithm is noticeably less expensive «han conventional 
error analysis methods. The U-D me.,;.od is shown to be especially attrac- 
tive for problems involving large numbers of bias parameters since 
it yields accurate and efficient techniques for performing sensitivity 
analysis and reduced-order filtering. 
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Chapter I . Introduction 


1.1 Background 

The optimal estimation of stochastic linear dynamic processes from 
imperfect linear measurements is a key problem in the fields of communica- 
tion and control. Usually estimators are linear functions of the data, 
and optimal solutions are those which minimize the mean square error. 

When estimates are based solely upon past and present measurements, 
this problem is termed a linear filtering problem. 

This report addresses the discrete -time linear filtering problem. 

A solution to this problem, for general nonstationary processes, was 
first derived by Kalman [I960]. The discrete Kalman filter is a recursive 
algorithm consisting of two parts: a time update and a measurement 

update. Each part contains difference equations for propagation of 
a state estimate and its error covariance matrix. The efficiency and 
simplicity of Kalman's algorithm make it particularly attractive for 
use in real-time estimation problems involving small digital computers. 
These features are among the reasons that Kali..\n filtering techniques 
have been widely used in a variety of engineering applications such 
as spacecraft navigation, aircraft guidance and control, marine navigation 
and power systems control (cf Battin and Levine [1970], Huddle [1969], 
Holdsworth and Stolz [1970] or Miller and Lewis [1971]). 

Although the discrete Kalman algorithm has been successfully 
employed in a number of filtering situations, practical applications 
of the method liave often been piqued with numerical difficulties. 


1 



Instances of serious accuracy loss in the Kalman filter have been reported 
by Bellantoni and Dodge [1967], Schmidt et al. [1968], and Dyer and 
McReynolds [1969]* Computational problems with Kalman's method are 
often evident in the form of Indefinite computed covariance matrices. 

Loss of covariance positivity is usually the result of computer round- 
off and cancellation' errors, aggravated by numerical ill-conditioning. 
Numerical problems may occur, for example, when very accurate measurements 
are processed in conjunction with large initial error covariances, 
or when a linear combination of parameters can be precisely estimated 
while others are relatively unobservable. In these cases computations 
involving the error covariance matrix are particularly susceptible 
to round-of and cancellation errors. 

A number of schemes have been devised to prevent loss of covariance 
positivity with the hope that algorithm performance would thereby improve. 
For example, Schmidt [1967 and 1968] has used artificially large process 
noise and measurement noise covariances, while Kaminski [1971a] suggests 
coordinate rotation. Such problem-dependent techniques are nonoptiraal, 
largely empirical, and often cumbersome. Moreover, they are usually 
inappropriate when precise parameter estimation is required. 

The computational shortcomings of Kalman's formula have motivated 
researchers to derive alternative formulations of his solution . While 


+ 

Significant cancellation errors can occur when two nearly equal numbers 
are differenced. Suppose, for example, that two numbers agree to 
six digits, and each number is accurate to eight digits. Then the 
difference would have only two-digit accuracy. 
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these methods are algebraically equivalent to the Kalman algorithm, 
they represent different computational techniques designed for improved 
numerical accuracy. The alternative algorithms usually require more 
c^putation than Kalman's original formula, and some methods Involve 
additional computer storage. Hence, many of these methods are not 
suitable for applications where computer time and storage are limited. 

This situation is common to many real-time applications, e.g., on-board 
navigation of aircraft or spacecraft. In the following paragiaphs 
several Kalman filtering techniques are discussed with this kind of 
application in mind. 

Alternative formulations of the Kalman algorithm generally fall 
into two categories: the information filters and the covariance filters. 

Information filters recursively compute either the information matrix 
or one of its square roots Covariance-related algorithms, like Kalman's 
original formula, deal directly with the error covariance or factoriza- 
tions of this matrix. Square root formulations in each category are 
acknowledged to be numerically superior to their conventional counterparts. 
Kaminski [1971bJ points out that square root factorization of a filter 
algorithm Improves numerical conditioning and provides greater effective 
precision. 

The square root Information filter was introduced by Oolub [1965] 
and Businger and Qolub [1965] as a reliable solution to the linear 


+ T 

The factorization As SS^ is not unique. See Appendix A. 
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least squares problem. This method, based upon Householder transforma- 
tion techniques (of Householder [1964]), was later extended by Dyer 
and McReynolds [1;-69J to include process noise. Their filter algorithm 
has been applied extensively in soacecraft navigation and has demon- 
strated superior numerical characteristics (cf Christensen [1976]). 

Analysis by Blerman [19733 and Kaminski [1971a] has shown that 
square root information filters are particularly efficient for problems 
involving large batches of data and infrequent estimate calculations. 
However, they show that covariance-type filters are more appropriate 
for real-time applications or when measurements are sparse. In these 
situations measurements are most efficiently processed one at a time 
(cf Bierman [1973]), in which case the covariance filter is referred 
to as a point processing algorithm. 

Recent research in covariance filtering methods has been stimulated 
by the need for fast, reliable point processing algorithms. A number 
of square root covariance methods have been investigated. The 
Chandrasekhar-type algorithms developed by Kailath [1974], Morf and 
Kailath [1975], and Lindquist [1974] appear to be efficient square 
root estimation schemes for stationary processes. These algorithms 
are not directly applicable, however, to the general nonstationary fil- 
tering problem considered in this report. For this reason Chandrasekhar 
methods are omitted from further discussion. 

Several covarianoe factorization algorithms have been derived 
to solve the nonstationary filtering problem. Notable among these 
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is the square root covariance filter, originally introduced by Potter 
[1963] and later extended by Schmidt [1970]. The Potter-Schmidt algorithm 
relies on Householder techniques which are known for their numerical 
stability and accuracy. Although the Potter-Schmidt filter has produced 
reliable results, Bierman [1973] has shown that this method requires 
considerably more storage and computation than does the original Kalman 
algorithm. This inefficiency is related to the fact that Potter's 
square root is a general nxn matrix, while the Kalman formula involves 
a symmetric matrix. 

Motivated by the need for a more efficient square root covariance 
filter, Carlson [1973] derived an algorithm which retains the square 
root in triangular form. Although sometimes less expensive than Potter's 
method, the n-dimenslonal Carlson filter requires n square root calcula- 
tions each time a scalar measurement is processed. Square root calcula- 
tions are usually time consuming compared to other arithmetic operations. 
Hence, for many applications Carlson's method is still noticeably more 
expensive than the Kalman formula. 

A promising new approach to Kalman filtering Involves a triangular 
covariance factorization which requires no square roots. The covariance 
P is factored such that P « UDU^, where 0 is unit upper triangular 
and D is diagonal. Bierman [1976a] suggested this factorization and 
derived a 0-D measurement update algorithm. The numerical integrity 
of this algorithm has been established by the work of Gentleman [1973] 
and [1975] which relates the U-D measurement update to the numerically 
stable Givens transformation methods. Moreover, the computational 
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requirements of the 0-D algorithm are equivalent to those of Kalman's 
method. Because of these attributes the U-D factorization algorithm 
is ideally suited for real-time applications. 

Bierman's results provided the starting point for this research. 
Attention was directed to the following problem areas: 

( 1 ) Extension of the U-D estimation method to allow for time 
propagation . 

(2) Demonstration, by analysis and experimentation, that the 
U-D factorization filter is a reliable and efficient point- 
processing algorithm. 

(3) Application of U-D filtering techniques to other areas 
of linear estimation theory. 

1.2 Outline of the Contents 

Chapters II and III contain a description of the various discrete- 
time covariance filter algorithms. In Chapter II attention is restricted 
to recursive data processing methods. The conventional Kalman measurement 
update formula is presented, and several alternatives to this method 
are described. These discussions address the computational aspects 
of each algorithm and conclude with a derivation of Bierman's U-D measure- 
ment update formula. Finally, a simple example problem is solved to 
lllustratec'the Improved performance obtained with the U-D and square 
root covariance factorization methods. 


6 


33-798 


Chapter III is devoted to time propagation. After the filtering 
problem is stated and the conventional Kalman solution is presented, 
several covariance propagation algorithms are described. These algorithms 
are designed to propagate covariance factors and generally Involve 
orthogonal transformations. Modifications of the familiar Gram-Schmidt 
and Givens triangularizatlon techniques are used to derive reliable 
U-D propagation algorithms. The general propagation methods are then 
adapted to problems Involving oias parameters and colored process noise. 
Exploitation of system structure yields a particularly efficient U-D 
colored npise time update algorithm. Finally, an example problem is 
Included to illustrate the superior numerical characteristics of the 
orthogonal transformation methods for covariance propagation. 

Chapter IV contains analytical cost comparisons of the various 
algorithms studied in Chapters II and III. Comparisons are based upon 
arithmetic operation counts which are weighted to reflect the different 
execution times required for each calculation. The measurement update 

p 

algorithms are compared first, followed by the time update cost compari- 
sons. Based upon this analysis the most efficient U-D and square root 
covariance propagation algorithms are selected. Measurement and time 
update costs for each method are then combined to yield filter algorithm 
cost comparisons. 

In Chapter V the U-D filtering method is extended to obtain flexible 
and concise algorithms for performing covariance error analysis. An 
efficient gain evaluation method, suitable for analyzing the effects 
of incorrect a priori statistics, is first derived. Analysis is then 
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extended to allow for evaluation of mlsmodeled data and state transition 
matrices. The general error analysis algorithm enjoys certain simplifica- 
tions when mlsmodeled bias and process noise parameters are evaluated, 
and these effects are described. The analysis of neglected bias parame- 
ters is also given special attention. In this case the U-D evaluation 
method allows for efficient sensitivity analysis and variable dimension 
filtering . 

In Chapter VI the various filter algorithms are applied to a 
realistic planetary navigation problem. Numerical accuracies of the 
different methods are compared by computing in both double and single 
precision arithmetic. Double precision results from all algorithms 
are in close agreement and are used as a reference for comparing the 
single precision results. Variations are introduced into the system 
model in order to evaluate algorithm sensitivity to a priori statistics 
and state dimensionality. Error analysis for this study is performed 
by applying the U-D gain evaluation method developed in Chapter V. 

Jc* 

Chapter VII gives a summary of results and recommendations for 
further research. 


8 




The notation x - N(x,P) describes x as normally distributed with mean x, 
and covariance P. Without loss of generality we assume the m components 
of v(i) to be uneorrelated . Correlated observations may be uncoupled by 
the "whitening* process suggested by Andrews [1968]. See Appendix A. 
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This thesis is primarily concerned with problems where m«n and 
where estimates are required frequently during the data processing. 
Bierman [1973] and Kaminski [1971a] have shown that this kind of estima- 
tion problem is most efficiently solved by recursive covariance-type 
algorithms such as the following Kalman formula. 

Let the variable denote the minimum variance estimate of x 

given data - [z(0), z(1),...,z(i-1)}. It has been proven^ that 

^(i-1) . where E is the expectation operator and "/" denotes 

the conditioning. Kalman [I960] derived the following formula for 
computing given and z(i). 

A(l) , A(i-1) ^K(l) (zu: - A(i)$<i-1)) (2.6) 

K(i) = P(l-1)A(i)^(A(i)P(i-1)A(i)'‘’ + R(i))"^ (2.7) 

^(i) s P(i-1) - K(i)A(i)t(i-1) (2.8) 


This recursion has the following initial values. 

x^’"'^ = y (2.9) 

P(-1) = P (2.10) 


^of Sage and Melsa [1971]. 
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The vector K(l) Is referred to as the Kalman gain and ^(i) is the error 
covariance matrix. Thus 


^(i) s - x)^} (2.11) 

Proof of this algorithm is given in numerous texts (cf Astrom [1970], 
or Sage and Melsa [1971]). We will not repeat these well-known proofs 
here but will, Instead, consider the computational aspects of Kalman's 
algorithm. Bierman [1973] has shown that the Kalman formula is more 
efficient when the m components of z(i) are processed one at a time. 

In this case the mxm matrix inversion in Bq. (2.7) reduces to a trivial 
calculation, and the cost of updating is a linear function of m. Since 
the measurement errors in each batch are uncorrelated , the data vector 
z(i) may be Included one component at a time by cycling through Eqs. 
(2.6)-(2.8) m times. This approach is emphasized and clarified by 
rewriting the Kalman formula as a scalar measurement update algorithm. 

For convenience we adopt the notation 


s . 

i . 

(2.1?) 

p . ?(i-i) 

t • ?(i) 

(2.13) 


The Kalman measurement update may then be written as follows. 
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Conventional Kalman Scal ar Measu.^ement Update 

^ = X + K(z - a^5?) (2.110 

K s f>a/(a^Pa + r) (2.15) 

^ = P - Ka^’p (2.16) 

Henceforth time dependence Is suppressed for notatlonal convenience 
unless required to avoid confusion. Where recursions are Involved we 
will rely on the superscripts and "a" to denote a priori and a 
posteriori quantities, respectively. 


2.2 Stabilized Kalman Algorithm 

Numerical difficulties with the conventional Kalman algorithm 
prompted Joseph to reformulate the covariance update (cf Bucy and Joseph 
[1968]). His method, referred to as the stabilized Kalman algorithm, 
computes ^ in the following way. 

^ s d-Ka'^) P d-Ka’’)’^ + KrK'f (2.17) 

This formula can be mechanized as follows. 

s I-Ka"^ (2. IP) 

«2 * “l ^ (2.19) 


IP 
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A second method, due to Bierman [19733, computes P using vector outer 
products in the following way. 


= Pa 

P^ = P-Kv^'^ 


( 2 . 21 ) 

( 2 . 22 ) 



V2 = P^a (2.23) 

P = (P^-V2K^) + KrK^ (2.24) 


The first mechod is not necessarily more reliable than the second 

A . , 

one even though P in Eq. (2. 20) appears to be a positive definite matrix. 
Furthermore, the original Joseph arrangement requires nearly an order of 
magnitude more computation than does the Bierman method.^ Because of 
this inefficiency and because there is no proof of improved stability, 
the first mech'.nization is omitted from further discussion. Even when 
the stabilized Kalman formula is implemented by Eqs. (2.21 )-(2.24) , it 
involves more than double the arithmetic operations required by the 
conventional Kalman method. 



'’’For the most part, these mechanizations involve addition and multipli- 
cation operations. The first method requires 1.5n^+2n^+n multiplica- 
tions, while the second involves only 4n^+4n such calculations. 
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Although generally more reliable uthan the conventional Kalman 
algorithm, the stabilize! Kalman formula is susceptible to numerical 
errors. This little publicized fact is demonstrated in Chapter VI. 


Potter [1963] observed that numerical problems associated with the 
Kalman algorithms were often evident in the form of Indefinite computed 
error covariances. In order to avoid such degradation he factored the 
covariance ^ so that 


A A/Vp 

P = SS^ 


(2.25) 


and derived an algorithm for recursively computing S instead of 
Potter noted that the covariance update for scalar measurements can be 
written as follows 


^ = f-Ka^P = §[I- - ff^lS*^ 

o 


(2.26) 


where 


f = S'^a 


(2.27) 


0 a r+f^f 


( 2 . 28 ) 


The initial factor, S^, may be uniquely defined by applying a Cholesky 
square root decomposition to P^ (see Appendix A). 
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By proper algebraic manipulation Potter found that if a constant y is 
chosen such that 


- X 

Y = X = 1/a (2.29) 

1 + yfi^ 


then 

[I-Xff’’] = [I-yXff’’] Cl-YXff'*’]'*’ (2.30) 


Equations (2.26) and (2.30) imply that S may be computed as follows. 


k = S-yKf'*' (2.31) 

K = XSf (Kalman gain) (2.32) 

The following efficient mechanization of the Potter algorithm 
was suggested by Bierman [1976a], 
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Potter Square Root Measurement Update 


= S^a 

(2.33) 

S 1/(r+f'^f) 

(2.34) 

= X/(1+VrX) 

(2.35) 

s 'Sr 

(2.36) 

= X + K[X(z-a^x)] 

(2.37) 

= ?-(YK)f’’ 

(2.38) 


This square root method guarantees positivity of the computed error 
covariance. It is also numerically better-conditioned than the Kalman 
algorithms. Bierman [1973] has shown that the Potter algorithm is 
equivalent to a particular Householder update (cf Appendix B). House- 
holder methods are known for their accuracy and stability, and so this 
equivalence establishes the numerical integrity of the Potter formula. 

Kaminski [1971a] suggests that when problems are ill-conditioned 
square root filters can provide twice the effective precision of covari- 
ance methods. In the case of the Potter algorithm, however, numerical 
stability is coupled with greater computational expense. Since the Potter 
update computes a general nxn matrix, his method requires nearly twice 
as much storage and calculation as the conventional Kalman update. 
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Because of this inefficiency the Potter algorithm is not attractive 
for use in many real-time applications. 


2.4 Carlson Triangular Square Root Algorithm 

Motivated by the stability and accuracy characteristics of square 
root estimation techniques and by the need for a fast, reliable point 
processing algorithm, Carlson [19733 derived an alternative square root 
formula. His method recursively computes an upper triangular covariance 
square root as follows.^ 


Carlson Square Root Measurement Update Algorithm 


f s s’^a f*'* = (f^,f2,...fn) 


(2.39) 


s r Kq s 0 


(2.40) 


For j = 1,2,...,n cycle through Eqs. (2.4l)-(2.45) . 


“j = ® j-1 + V 


(2.41) 




(2.42) 


"'J = 


^The algorithm given here is a modest rearrangement of C arlson' s formula. 
Bierman [ 1976b] has observed that the calculation of in 

E q. (2.42 ) may be more accurate than the corresponding computation, 

V ®J-1 ’ '*®coroo®*'8ed by Carlson. 
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* ®J^j * ^j^J-1 


(2.44) 

Kj - ■^J-l ^ f J§j 


(2.45) 

vrtiere 



S — and — ^^J ^ ^ ^ * 

.,Sj(J),0,...,0)'”’ 


The Kalman gain is given by 




K = K„/on (2.46) 

Thus, the estimate update may be written as follows. 

X = X +lKp[(2 - a%)/otn] (2.47) 

Derivation of this algorithm is deferred until the next section 
where it is shown to be an easy consequence of the U-D measurement 
update . 


The Carlson algorithm has a good computational form and enjoys the 
traits of stability and accuracy generally attributed to square root 
filters. Although the Carlson formula requires considerably less storage 
and computation than the Potter square root method, It is still noticeably 
less efficient than the conventional Kalman algorithm. Unlike the conven- 
tional method, the Carlson algorithm requires n square root calculations 
for each scalar measurement update, and square roots are usually 
time-consuming operations. This is particularly true for small on-board 
computers such as the Litton 4516. This computer has 32 bit double 


BISBOPUCIMIjOT 0^ 
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precision arithmetic and requires approximately 1000 usee to calculate 
a square root. By contrast it performs an addition in only 4 usee. 

A careful comparison of measurement update costs is included in Chapter 
IV where the Carlson method is shown to be unreasonably expensive for 
a significant class of problems. 

2.5 U-i-D Covariance Factorization Algorithm 

Bierman [1976a] recognized that square root calculations required 
by the Carlson algorithm are often costly and proposed a square-root-free 
measurement update scheme. His method employs the covariance 
factorization 


P s UDU^ (2.48) 

where U is unit, upper triangular and D is a positive diagonal. It is 
well known that for symmetric positive definite matrices thia factoriza- 
tion exists and is unique (cf Martin et al. [1965]). An algorithm 
for computing the U-D factors of P is described in Appendix A. 

Measurement updating using the U-D factorization preserves the 
non-negative definite structure of P and is equivalent to the Carlson 
method without square root computations. 


^The factor U could also be taken as lower triangular. However, upper 
triangular factorizations allow for variable dimensioned filtering 
as described in Chapter V. 
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The U-D data processing slgorlthm may be derived by factoring the 

^ A 

Kalman update, Eq. (2.16). IfP and P are both factored according to 
Eq. (2.48), then 


fi & = U[D - - w'‘’]u'*‘ 


(2.49) 


where 


V = Df = djf 


(2.50) 


f = U^a 


(2.51) 


+ fT^f = r + £ v^f 
1=1 


The bracketed term In Eq. (2.49) Is positive definite and hence may 
be factored as U D U^. Since the product of two unit upper triangular 
matrices Is again unit upper triangular it follows that 


U = U 0 6 s D 


(2.53) 


Hence the U-D update rests on the I'actorization 


U 6 0"^ = D - (1/«) vv**^ 


(2.54) 


Bierman has observed that the Agee-Turner matrix factorization in 
Appendix C may be applied to this problem and results in the following 
recursion for computing 0 and 6. 
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With Cp s -1/«, cycle through Eqs. (2.55)-(2.57) for j = n,...,2. 


^3 2 

dj . dj ♦ Cjv/ 


“kj = ’ 




d^ s d^ + c^v<j 


Equation (2.55) is subject to cancfcllation errors since it 
involves the differencing of two positive quantities. However this 
expression for Sj- and, hence the entire recursion, can be rewritten 
as a numerically stable formula. 


Let the partial sums, Cj, be defined as 


®j = ■' + i ^i^i^ = ®J-1 


Then since Vj * ^j^J» ®9®* ^2.55) and (2.57) imply that 


Oj . -l/.j 


This value for Oj yields the following expression for dj. 


■ ■. {¥) 


J * 1,2,...,n 


(2.61) 
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Sqs. (2.56}-(2.6l) imply that 

IT s I ♦ [0lX2V^^^IX3V<2)j _j (2.62) 

trtiere 

. (v,,Vj Vj, 0,0,0) (2. 64) 

The desired factor 0 s U IJ may now be constructed from Eq. (2.62). 
Let Oj and Uj denote the Jth columns of U and U respectively. 

Then 

= ®J * ‘A-1 <2-«5) 

Tj = 0.'J> = ICj., . VjOj (2.66) 

'K'o = 0 (2.67) 

and 

Xp/o^j = K (Kalman gain) (2.68) 

The U-D update Is summarized as follows. Given a priori covariance 
factors U and S and the scalar measurement zsa^x>v (E{v^)sr), the 
updated covariance factors 0 and B and the Kalman gain, K, are obtained 
by evaluating the following sequence of equations. 
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Blerman U-D Measurement Update Algorithm 


f s U^a fT = (fpf2,...,fn) (2.69) 

V s Df Vj = dj fj (2.70) 

n-1 

k/ s (v^,0,...,0) (2.71) 

s r + (2.72) 

s d, (2.73) 

For J = 2,...,n cycle through Eqs. (2.74)-(2.78) . 

*j = °j-i * ’/j 

• (^) 3j <2-75) 

‘j . (2.76) 

. Oj . ‘/j., (2.77) 

*J • *J-t * ’J^J 

where 


0 * [Ui,O2.-*-0n3 and 0 « [0, .Og,. ..8n] 
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The Kalman gain Is given by 



(2.79) 


The salient feature of Bierman's algorithm is the way in which the 
updated diagonal B is computed. Since the quantities, aj, are calculated 
as positive sums (cf Eq. 2.74), cancellation-type errors are avoided and 

A A 

the positivitj of D, and hence, P, is assured. Furthermore, the elements 
of D may diminish to near-zero without affecting the stability of the 
algorithm. The numerical integrity of this method is further established 
by the analysis of Gentleman [1973] and Fletcher and Powell [1974]. 

Their work relates the U-D algorithm to a numerically stable Givens 
orthogonal transformation. An efficient Fortran implementation of 
the U-D measurement update algorithm is included in Appendix D. 

Proof of the Carlson algorithm (Section 2.4) is immediate from the 
U-D formula when the square roots S and S are identified as 

5 . Oo’'^ § = (2.80) 

Equation (2.80) suggests that the U-D method shares the attributes 
of accuracy and stability generally associated with square root filtering 
techniques. This fact is demonstrated by the case study in Chapter VI. 

2.6 

The following example provided by Blerman [1976b], illustrates the 
numerical characteristics of the various data processing techniques 
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described in this chapter. In particular it demonstrates how the conven- 
tional Kalman method can produce nonpositive error covariances while the 
factorization methods yield accurate results. 


Consider the problem of estimating and %2 from scalar measure- 
ments z^ and Z 2 vdiere 


z^ s ♦ ex2 + v-| 


Z2 5 ♦ X2 ♦ V2 


( 2 . 81 ) 


Data errors, v.| and V 2 are uncorrelated, zero mean random variables with 
unit variances. The a priori error covariance is assumed to be 


P S 0^1 

where os1/e and 0<c«1. The quantity e la assumed to be small enough 
such that computer round-off produces 


1 ♦ c2 * 1 


( 2 . 82 ) 


This estimation problem is certainly well posed. The observation 
z^ provides an accurate measurement of x^ which, when coupled with the 
observation % 2 t should accurately determine X 2 < However, when the 
various data processing algorithms are applied to this problem several 
diverse results are obtained. 
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Let the gain and error covariance associated with the measurement 
Zf be denoted as K^ and P^ respectively. Similarly the measurement Z2 Is 
associated with K 2 and p2< Table 2.1 gives the exact solutions for K and 
P at each step and shows the rounded solutions computed by the various 
filter mechanizations. Notice how computer round-off errors are evident 
after the first measurement is processed. The first covariance matrix 
computed by the conventional Kalman formula has a zero variance. One 

p 

might expect this result if the exact variance were on the order of c or 
even e. Howtsver, in this case the correct value is approximately = 2, 
and so the Kalman error is appreciable. Subsequent processing of data 
by the conventional formula results in even larger errors. The second 
computed covariance, for example, has one negative diagonal element 
and one that is zero. It is also interesting to note that this covariance 
matrix has off-diagonal elements equal to -1 if the P12 element is 
computed and set equal to P2^. On the other hand, if P21 is computed 
directly it has a value of 4I. Yet computer Implementations of this 
algorithm typically Include calculation of only the upper (or lower) 
triangular elements. 


Although the stabilized Kalman formula performs better than the 
conventional method, it too computes different values 'or P^2 ^21 

after the second update. The quantities P12 - 3c and ?2^ s -1 - 2c 
cculd be averaged to give P^2 ~ ^21 = " (5/2) e. However, note that 

the exact answer rounds to P^2 ~ ^21 ~ ~ 3c. 
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All the factorization algorithms discussed here produce accurate 
results. The U-D and Carlson methods yield nearly identical results 
and compute rounded covariances equal to the rounded exact answer. 

This sample problem illustrates how the Kalman covariance algorithms 
a*’e prone to numerical errors which can significantly affect filtering 
accuracy. The factorization methods, on the other hand, avoid critical 
round-off and cancellation errors and yield reliable results. 
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Table 2.1. Comparison of Solutions to Example Problem 



•Covariances are not a part of the Potter, Carlson and U-D algorithms. 
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Chapter III . Filtering for Discrete Linear Processes 

3.1 Problem Statement and Conventional Solution 

The discrete linear filterliig problem is an extension of the 
parameter estimation problem discussed in Chapter II. Instead of being 
constant, the parameter x is described as the state of the following 
linear multistage process. 

x(i+l) = •(i)x(i) + 3(i)w(i) i = 0,1,2,... (3-1) 


I 


i 



I 



I 








The array '*(1) is a known nxn transition matrix, B(i) is a known 
nxk matrix and 


x(0) - N(x,P) 


(3.2) 

w(i) - N|o,Q(i)j 

Q(i) = diag (Qi , . . , q^) 

(3.3)‘'‘ 

E{w(i)w^(j)} = 0 

i j 

(3.H) 

E{x(0)w^(i)) = 0 

for all i 

(3.5) 


At stage 1 scalar measurements are available and are linearly related to 
x(i) by the equation 


^There is no loss of generality in assuming that Q(i) is diagonal. 
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z(i) = a^(i)x(i) + v(i) (3.6) 

where 

v(i) - N(o,r(i)j (3.7) 

S(v(i)v’‘(j)} =0 i ^ j (3.8) 

E{v(i)x^(j)} = 0 for all i,J (3.9) 

E{v(i)w'‘^U)) = 0 for all i* j (3.10) 


The optimum solution to this filtering problem is the minimum 
variance estimate of x(i) given the observations = {z(0/ ,z( 1 ) , . . . ,z(i)} . 
Kalman [I960] derived a recursive solution to this problem by combining 
the data processing formula of Chapter II with an optimal time update 
algorithm. The complete Kalman filter is summarized as follows. 

Conventional Kalman Filter 

♦(i) ^(i) (3.11) 

♦(i) P(i) ♦'''(i) + B(i) Q(i) B^'U) (3.12) 


| x(i+1) = 
P(i+1) = 
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x(i) = x(i) + K(i)(z(i) - a’’(i) x(i)) 


(3.13) 


Measurement 

Update 


^(i) = P(i) - KCDaTcDPU) 


(3. 1^0 


K(i) = P(i)a(i)/(a'*‘(i)P(i)a(i) + r(i)) (3.15) 


This recursion is initialized by the following quantities. 


x(0) = X 


P(0) = P 


( 3 . 16 ) 


In Chapter II several alternative formulations of Kalman's measure- 
ment update algorithm were presented. These data processing algorithms 
may be combined with appropriate time update schemes to provide alterna- 
tive mechanizations of the Kalman filter. For example, the stabilized 
Kalman measurement update, Eq. (2.17), may be substituted for the conven- 
tional Kalman formula, Eq. (3.111), to yield the stabilized Kalman filter. 


Measurement update algorithms involving covariance factorizations, 
such as the Potter square root or the Bierman U-D methods, require alter- 
native formulations of the covariance propagation formula, Eq. (3.12). 
Several techniques for propagating covariance factors are described 
in the remainder of this chapter. 
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3.2 RSS Method for Propagating Covariance Factors 


Suppose that at each stage of the filtering process the a posteriori 


covariance, P, is factored such that 


^ or ^ = UD 


(3.17) 


The RSS method for computing the a priori factors S (or U-D) at the next 
stage involves "squaring up" S (or 0-D) to obtain The conventional 
propagation formula Eq. (3.12) is then applied to compute P, followed by 
an appropriate Cholesky decomposition to obtain S (or U-D). 


Carlson [19733 recommends the RSS method as an efficient square 
root propagation scheme. However, the analysis in Chapter IV demon- 
strates that this formula enjoys only a minimal savings in computation. 
Of greater significance is the potential loss of accuracy associated 
with this method. A major motivation for factoring the Kalman algorithm 
is to gain increased accuracy due to better numerical conditioning. 

This advantage can be eliminated by applying the RSS propagation scheme. 
The example problem in Section 3.7 Illustrates why the RSS method is 
numerically hazardous. 


3.3 Gram-Sohmidt Propagation Algorithms 


Suppose ^ is factored as P = and let 


W = [B I *6] f n 


( 3 . 18 ) 
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k n 

D = diag (Q, 5) (3.19) 

The propagation equation P s + BQB^ may then be written as 

P = WOW*** (3.20) 

The U-S factors of P may be obtained by applying a particular Gram 
Schmidt orthogonalization to the row vectors of W. This procedure 
involves the following inner product space. 

D 

Let Vjj denote the vector space of n-tuples over the reals with 
inner product, <.,.>d where 

fp D 

<Vi,Vj>D = Vj^^DVj, for. all v in 
||v||p = <V,V>p 

Thus, the vectors v and w are "D-orthogonal" if <v,w>p = 0. 

n u 

Suppose is an independent set of vectors in V^j where n i m. 

A D-orthogonal set {vj^} may be obtained from {Wj^} by applying the follow- 
ing weighted Gram-Schmldt (WGS) algorithm. 


( 3 . 21 ) 

( 3 . 22 ) 
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Mfelghted Gram-Schmidt Ortfaogon alizatlon 


’'n-1 = «n-1 


Vd 


ii-nii; 


2 n 


vi = w, 


E -ll 

I Iv 


<Wi,V.> 


i'D 


i=2 


2 ’i 

i"D 


IvJI 


(3.ir3) 


This WGS algorithm may be used to compute the U-D factors of 

P in the following way. Consider the n independent row vectors of 

D - 

W in Eq, (3.18) as vectors^^ where m = n k. If the D-orthogonal 


set is computed via applied to then Eqs. (3.23) imply 



^ (3.2i|) 


where U is unit upper triangular. Since sectors {vj^} are D-orthogonal 
it follows that 
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P = WDw"^ = U 

Thus, the elements of U and D are given by 

2 

dj = llVjllp j = 1,...,n (3.26) 

1 

Uij = — <Wi*Vj>5 i = (3.27) 

The orthogonalization used to compute U and D Is a modest generali- 
zation of the familiar Gram-Schmidt (GS) procedure (of Noble [19693) 
in that inner products are weighted by the D matrix and the vectors 
{Vj^l are not normalizev. . The required orthogonalization may also be 
attained by applying a Modified Gram-Sohmidt (MGS) procedure, (cf Lawson 
and Hanson [197*1]). BJorck [1967], Jordan [1968] and Rice [1966] have 
investigated the numerical characteristics of GS and MGS. Their studies 
establish that MGS is more accurate, takes no more arithmetic operations 
and requires less storage than does GS. Moreover, the works of Bjorck 
and Jordan show that triangularization using MGS has accuracy that 
is comparable with the reliable Givens and Householder methods to be 
described in this chapter. For these reasons the following Modified 
Weighted Gram-Schmldt (MWGS) algorithm is recommended for U-D time 
updating. 


Ivillp 


living 


UT 


(3.25) 


' Vd 
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Modified Weighted Gram-Sehmldt Factorization 

Given the full rank nxm matrix W with row vectors {w^} and an 

m-dimensional positive diagonal matrix, The U-D factors of 
~ — T 

P = WDW are computed as follows. 


Evaluate Eqs. (3.28)-(3.30) recursively for j = n, n - 


I ^ 


C 


dj = Ilwj<"-J)|l5 


(3.28) 




^ij = r- < > 5 




i = 


(3.29) 


(3.30) 


^ 4 


This recursion begins with the vectors where 


» ^ 


= Wj^ i s 1 , . . . ,n 


This algorithm may be obtained from the WGS results with the 
easily proven identity 


I .. 

* 'f 

} t 


Vi = 


'J “ "J 


J * l,...,n 


(3.31) 


Although the inclusion of superscripts makes the MWGS algorithm appear 
rather complicated, this method is easily mechanized. Compactness, 
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efficiency and simplicity can be achieved by arranging the FORTRAN 
implementation so that successive w's in Bq. (3.30) are written over one 
another (see Appendix D). 

Gram~Schmidt orthogonalization may also be used to propagate 
covariance square roots. Let P be given by Eq. (3.12) where ^ 

A triangular factorization P = SS^ may be obtained by the following 
MGS algorithm. 

Modified Gram-Schmidt Square Root Fact orization 


Let 

k n 

W = j ^3| n (3.32) 

with row vectors denoted as An upper triangular factor, 

§, such that 


(3.33) 

may be obtained from the following recursion with inner products 
defined by 


<Wi,Wj> = Wi^Wj 

Fcr J s n,...,1 cycle through Eqs. (3.35)-(3.39) . 


(3.34) 
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(3.35) 




(3.36) 


vj = 


JJ 


(3.37) 


®ij = vj> J (3.38) 

I i = 1,2 j-1 

y^(n“J+l) s - Sj^j Vj) (3.39) 

This algorithm is an easy consequence of the NWGS factorization, 
Eqs. (3.28)-(3. 30) , and the identities 


D = I S = 


(3.40) 


The MGS formula may be used in conjunction with both the Potter and the 
Carlson measurement updating methods. 

The recursions Eqs. (3.35)-(3.39) may be used to compute a lower 
triangular S if the indices, i and J, are reordered so that J * 1,...,n 
and i s J,...,n. Since the Potter algorithm does not require § to 
be upper triangular, this rearrangement of the MGS formula is also 
appropriate for propagating the Potter square root. 
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3.0 ClYens Transfonnatton Methods for Covariance Propagation 

Consider the problem of constructing an mxm orthogonal transforma- 
tion T such that 

m-n n 

[W]T s n jcTlIir] (3.41) 

where V is upper triangular. If W is defined to be 

k n 

W = n| I (3.42) 

then the transformation in Eq. (3.41) represents a covariance square root 
time update with V' s s. 

Givens [1959] showed that T could b.? constructed as a product of n 
orthogonal transformations, Tj, where 

T = T„ Tn_i...Ti (3.43) 

Bach Tj is designed to zero out the subdiagonal elements of row J. 

Thus, Tjj is constructed such that for a w we have 
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m-1 

1 




n {[W<"hT„ = = n-lj 


• 

• 

• 

JL 


or 0 0 o' 

• 

— 



(3.410 + 


This transformation consists of a sequence of two~dimenslonal column 
rotations called Givens transformations. These rotations pivot in 
succession the last column of the array H with each of the preceding 
columns. The rotation Involving column 1 is designs .1 to set * 0- 


In a similar manner is consti*ucted to zero out the subdiagonal 
elements of row n-1 by successively pivoting column m-1 of with 

each of the preceding columns. Notice that this transformation does 
not alter the last row and column of This process is continued, 

and after Tj^^ is applied the array has the form 



(3.l^5^ 


^An asterisk is used to denote 


nonzero, unnamed elements. 
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Thu3, the process is concluded with the transformation T^, 

This triangularization technique involves a sequence of elementary 
Givens rotations of the following form. 


yi 

^2 


n-j 


0 


T = 



(3.46) 


% ;; 
I " 
I 1 

I ; 

I ^ 

V 

I. 

I' = 


The vector y is the pivotal column and 


yj = * yj^ 

(3.47) 

° = yj/yj 

(3.48) 

s = Xj//3 

(3.49) 

= cxj^ - syj^ ' 

(3.50) 

( i=1....,J-1 


y^ = sxj + cyj ) 

(3.51) 
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The Givens covariance square root propagation algorithm is suinnarized 
below. Superscripts have been suppressed to simplify the notation and 
to facilitate computer implementation. The symbol is used instead 
and denotes replacement in computer storage. 

Givens Square Root Fact orization Algorithm 

Let V be a full rank nx(n+k) matrl- with column vectors 
{Wi)i_^. An upper triangular factor S such that SS^ =ww'*’ may be computed 
as follows. 

For J = n,...,1 cycle through Eqs. ^3.52)-(3.58) . 

m := k+j (3.52) 

For i = m-1,...,1 evaluate recursively Eqs. (3.53)-(3.58) . 



( 3 . 53 ) 

c := W„(j)/w'(j) 

( 3 . 54 ) 

s := Wi(j)/w^(j) 

( 3 . 55 ) 

■H 

II 

> 

( 3 . 56 ) 

Wi := cwj_ - sw„ 

( 3 . 57 ) 


k2 
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“m (3.58) 

k n 

When this recursion is completed the W array has the form W = [0 ! S][n 
where S is upper triangular. 


A lower triangular factor S may be obtained from the recursion 


Eqs. (3-52)-(3.58) if the indice^T^and j, are reordered so that 
j = 1,...n and i = j+1,...,m. At the conclusion of these calculations 
the W array has the form 


W = [S ! 0] f n 
where S is lower triangular. 


(3.59) 


The Givens triangularization method is well known for its superior 
numerical characteristics (cf Wilkinson [1965]). Thus, the Givens 
square root factorization algorithm provides i reliable method for 
propagating Potter and Carlson covariance factors. However, this method 
is usually bypassed in favor of the more efficient Gram-Schmidt and 
Householder time update algorithms. The reader is referred to Chapter 
IV for detailed cost comparisons. 

A modest generalization of the Givens triangularizat ion technique 
yields a reliable and efficient U-I) time update scheme. Recall that U-D 
factorizations involve the propagation equation P = WDW^ where 


1^3 
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k n k n 

W = [B I *U] } n and D = diag (Q,^ 

Let {Wj^} denote the column vectors of W and let {dj^} represent the 
diagonals of D. Note that P can be written as 

( 3 . 60 ) 

^ Wg I ... ! Wjij , msn+k (3.61) 

The Givens triangularization method may be applied to the array W in a 
way which explicitly accounts for column scaling. That is, an orthogonal 
transformation T may be constructed such that 


P = HW^ 

where 



m m-n n m-n n 
n {[W]T = [ 0 ! ^ u^ ! ... I u„] = [ 0 I (3.62) 

where U is unit upper triangular. Genbleman [1973] showed that an 
elementary Givens rotation, adapted for such column scaling, takes 
the following form. 


Itl* 
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/axy /ff yJ 


^'x\ /^y 


where 


/ByJl 


2 o 2 

e xj + B yj 


= Bo/B^ 


c S By^/B"' 


s = <Wj/8^ 


*i = yj*i - *jyi 


y^ = oy^ + sx^ 


y.i-1 


i=1,...,J-1 


These expressions for o", B"*, x' and y' are easily obtained by applying 
Eqs. (3.47)-(3.51 ) with x and y replaced by /a x and /ff y respectively. 


Note that when the quantities a , B , x and y' are desired, 
rather than the vectors (4^x0 and this method involves no 
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square root calculations. Gentleman [1973] noted that when yj = 1, 
Eqs. (3.6*0-(3.69) can be rewritten to avoid a multiplication. Thus 


X 2 

B B ♦ «Xj (3.70) 

a' = fla/e" (3.71) 

s = axj/B' (3.72) 

= *i - XjYj (3.73) 

y£ = 7i + s:t£ (3.74) 


This algorithm is referred to as Method B while the formula described by 
Eqs. (3.64)-(3.69) is denoted as Method A. Both of these methods may be 
used to solve the U-D time update problem. 

The U-D factorization of P = WDW^ is initiated by denoting as 
pivotal elements the pair (d^iWj,) where w^ represents the last column 
of W. The square-root- free Givens transformation defined as Method 
A is then applied to the pairs (dj,_p Wn-1^ ^^m» '‘m^* ® result 

of this initial pivot w^_^(n) = 0 and w'(n) = 1. Thus subsequent 
transformations Involving the mth column may be accomplished by applying 
Method B. The pair (d^,w^) is then pivoted in turn with each of the 
remaining pairs (dj^,Wj^', 1 < m-1. Implicit in this sequence of calcula- 
tions is the transformation Tjj where 


46 




= diag (d(, d 2 ',...,d^.i, (3.77) 

This technique is repeated, with the pair (d'_^,WjJ_^) as the pivotal 
elements, to yield and The process is thus continued 
until finally and are obtained where 


m-n n 
= [ 0 I U] } n 


1»7 


(3.78) 





f 
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I 

4 


m-n 


= diag (•••», d^,d 2 , 


(3.79) 


The modified Givens algorithm for propagating U-D covariance factors 
is sununarized below. 


Modified Givens Factorization 

Suppose we are given the full rank nx(n+k) matrix W with column 
vectors {Wj^} and an (n-t-k)-dimenslonal positive diagonal matrix, D. 

The 0-D factors of P = HDW^ may be computed as follows. 

For J s n,...,1 cycle through Eqs. (3.80)-(3.89) . 


m := J+k 


( 3 . 80 ) 


For i s m-1,...,1 evaluate Eqs. (3.8l)-(3.89) as indicated. 

:= + ^i (3.81)^ 

s := w^(J)/d^ (3.82) 

^i •= W'lm <3.83) 


1 

I 




s' 


I 

i . 

I 

'I 

i 

i 


■^When i < m-1, Wj^(j) s 1. 



& 


f' 


1*8 
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'= ^*111 if i < m-1 


(3.810 


(3.85) 


”i •= "i " Wi<J)Wn 


( 3 . 86 ) 


Method A: Wj^ ;= cWju + 5v if i = m-1 


(3.87) 


Method B: := Wj^ + swj^ if i < m-1 


( 3 . 88 ) 




(3.89) 


Upon completion of this recursion the W array contains U. That is 


k n 


W = [0 ! U] }n 


(3.90) 


The diagonal elements of D are given by 


d^ - djj^^ i s 1,,..,n 


(3.91) 


Analysis by Gentleman [19733 has established that when Method A, 

Eq. (3.87), is used exclusively in this recursion the algorithm is always 
numerically stable. However, use of Method B, as indicated in Eq. (3.88), 
saves approximately 1/3 n^ multiplications (33*). Hence this formula is 
preferred when it can be relied upon. Gentleman [1973] has shown that 


1*9 
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errors associated with Method B can be related to the ratio d;;/d_. 

He recommends control of error growth by testing this ratio against some 
limiting value, c. Thus, Method B is applied only when d'/d^ < c. 
Fletcher and Powell [197^] have experimented with this technique on a 
large number of problems using c = 4.^ Their experience suggests that 
for most applications of this algorithm Method A is required so rarely 
that efficiency is near optimum. A numerically reliable computer luple- 
mentation of the modified Givens update is Included in Appendix D. 






Recall that propagation of covariance square roots involves the 
triangularization problem. 


m-n n 


n { [W] T = n I [ 0 IS] 


where S is upper triangular and T is an roxm orthogonal transformation. 
This problem may be solved by applying a sequence of Householder trans- 
formations, Tj, such that T = TjjTjj_^ . . .T^ (cf Householder [1964]). Each 
elementary transformation is a reflection operator and may be written as 


vrtiere 


Tj s I - Rj 


Bj = 2 /||u^J^||2 


(3.92) 


(3.93) 


tGentleman suggests c s 10 as a reasonable limit 
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Note that for any vector y 
y^Tj = 

with Y given by 
Y = 


(3.94) 


(3.95) 


Thus, explicit formation of the matrix Tj is not required if only the 
product Tj is desired. 

The triangularization of W is begun by designing the transformation 
Tj^ to zero out the subdiagonal elements of row n. Thus, 


[W]Tn = = 


n-1 


I 


m-1 

1 



w' 

^2 


'^r-1 

0 0 0 0 

"°n 




(3.96) 


Equations (3.92)-(3.96) imply that T„ can be constructed from the vector 

u(h) . jj(n) ^ 


«n«m 


(3.97) 
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where 


is the nth row of the array W and 


m-1 

s (0,0, 0,0,1) (3.98] 

Since Tjj is orthogonal, is constrained such that s ||x^”^||^. 
Thus Eqs. (3.93) and (3.97) imply that 

»n = '/(•/ * <3-«: 

In order to avoid cancellation errors when is calculated we may 
define 


= sgn (x^"^)||x^"^|| (3.10( 

Once Tjj has been applied *^o the array W via Eq. (3.9*0 » the trans- 
formation is similarly chosen to pack the (n-l)st row of the array 
y(n-1)^ Construction of is identical to the construction of u^” 

with one exception. The mth component of is free to be chosen 

since there are only m-2 subdlagonal elements to be zeroed out. If the 
mth component of is set to zero then (0 0...0,-Ojj) « 0. 

Thus, Eqs. (3. 9*>)-(3.95) imply that T^^. will not disturb the nth row 
or the mth column of the array. Hence, 


RIPRODUCIBILITy OP THi 
18 POOB 




n-1 


( 3 . 101 ) 



This process is continued so that at step j is constructed 
according to Eq. (3.97) with n replaced by J and m replaced by (m-n+J). 
The m-dlmensional vector is defined as follows. For j = n,...,1 



1 = 1 , . . . , J+m-n 


i > J+m-n 


where w 


(J) 

J 


is the Jth row vector of 


( 3 . 102 ) 


The following algorithm, due originally to Schmidt [19703, applies 
the Householder triangularization method to propagate covariance square 
roots. 


aohMldt COYtrianoe Scuare Hoot Time nndate 

Let P ■ and define W to be the nx(n+k) array 
k n 


W « I jn 




with elements An upper triangular factorization, SS^ s ww", 

may be obtained as follows. 

For j » n, n-1,...,1 compute recursively Eqs. (3*103)-l3.111). 



6j * 


For isl ,2, . . . , j-1 cycle through 

J+k 

n •’ «it 

Isl 

«ii •* «il - n 

•* -®J 

«Jl •* 0 


sgn (Wj,j+k^ (3.103) 

i = (3.104) 
i = J+k (3.105) 
i s j+k+1,...,n+k (3.106) 


(3.107) 

Eqs. (3.''C8)-(3. 109). 

( 3 . 108 ) 


ts1 J+k (3.109) 

( 3 . 110 ) 

i s i,...j+k-l (3.111) 




im 
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I 

I At the conclusion of this recursion the array W has the form 

I k n 

I W = [0 ! S]| n 

where S is upper triangular. 

The Householder trlangularlzation method can be modified to yield 
a U-D propagation algorithm by applying the identity S = OD' to 
Eos. (3.103)-(3. Ill ) . However, this method is not pursued sin- e it 
appears to have little advantage over the MHGS and modified Givens 
algorithms. 

3.6 Propagation Algorithms for Systems with Colored Process Noise 

Time propagation represents a major filtering expense, particularly 
in applications where measurements are sparse (cf Chapter IV). Hence, 
time updating methods warrant further consideration in terms of efficient 
computer implementation. It is often possible to significantly reduce 
propagation costs by exploiting special system structure. For example, 
one may trim computations considerably by taking advantage of sparse 
state transition matrices or triangular system equations. This section 
Illustrates how triangular covariance factorizations lend themselves 
to efficient propagation schemes for systems with bias parameters 
and colored process noise. Suppose ...e state transition is given by 


V 

r* 

X 

♦ 1 
xp 


X 

+ 

o' 

p 

J 

M 


p 


w 


i+1 L 

. 


. . 

i 

.. • 


i=0,1,... (3.112) 
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M — t * * * 


(3.113) 


E{w^ Wj } = Q Q = diag (qi,...,q|() 


ECxqWj'*'} = 0 E{PqWj^'^} = 0 


This system model closely approximates a large class of second- 
order processes and hence has wide application (cf Christensen [1976] 
or Maughmer and Byrd [1969]). The structure of Eq. (3.112) permits 
time propagation to be performed in two phases. Thus, 


*x *xp1 p 


I P 


0 1 rx 


M p 


Equation (3>115) defines a deterministic pha&e of the mapping. 
The error covariance associated with this phase is denoted as P, while 
the final map, Eq. (3.116) yields the a priori error covariance, P. 


Tor notational brevity we use subscripts in Eqs. (3.112) - (3.116) 
to denote time dependence of system parameters. Notice further that 
the matrices *jjp, M ‘».nd Q may also be time-varying. 
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The original problem is decomposed in this manner because each subproblem 
enjoys certain computational advantages, particularly when upper triangular 
covariance factorizations are involved. 


3.6.1 D-D Colored Noise Time Update 

Consider the covariance factorization P = UDU^. The simplifica- 
tions associated with a deterministic map of U and D via Eq. (3.115) 
are illustrated in the following algorithm. 

Bias Partitioned 0-D Factorization 
T 

Let P = UDU vrtiere 0 and D are given by 


U = 


xp 




D = diag (D^, Dp) 


(3.117) 


Assume P is given by 


P = *P*^ 


( 3 . 118 ) 


where 


» s 


'xp 

I 


I" 


(3.119) 
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Then the U-D factors of P, with similar partitioning assumed, are deter- 
mined as follows. 


[Op. Op] = [Op, Bp] (3.120) 

^xp ~ ^xp *xp ^p (3.121) 

VxV ■ >W B, <3- '22) 

Equations (3. 120)-(3. 122) may be obtained by substituting the 
partitioned U-D and ♦ factors into Eq. (3.118) and expanding the result. 
Note that Ujj and moy be computed by applying either the MWGS algorithm, 
Eqs. (3.28)-(3.30) , or the modified Givens algorithm, Eqs. (3.80)-(3.89), 
to the arrays W = and 5 = 0^^. 

An easy generalization of tnis partition algorithm involv replace- 
ment of the identity matrix in Eq. (3.119) with a nonsingular, t 
triangular matrix, *p. The modifications required for Eqs. (3.120) 
and (3.121) are apparent. 

The following features of the bias partition algorithm are notable. 

1) The U-D factors corresponding to the bias portion of the 
state vector remain constant. This property is contingent 
upon having the bias parameters in the lower portion of 
the state vector. 


58 


Rli3»R0DUCIBILITY OP THE 
ORKHNAL PAGE IS POOR 





33-798 

2) The cross coupling U„_ is updated by a simple linear relation. 

xp 

3) The updated factors 0^^ and are obtained by considering 
only that portion of the problem which is independent of p. 

4) All of the above comments also apply to problems where 
X is a stochastic process. In this case the right hand 

side of Eq. (3.122) would have an additional term correspond- 
ing to the process noise disturbance. 


Once the deterministic time update, Eq. (3.115), is accomplished 
with updated covariance factors, U and D, computed from Eqs. (3.120)- 
(3.122), the second phase of the mapping, Eq. (3.116), is performed. 

The diagonal structure of M and Q may be exploited by mapping the process 
noise, p, one component-at-a-time . The following algorithm indicates 
how the intermediate factors, U and S, are updated as each component 
of p is mapped. 

Single ComDor.ent_U- D Time Update 


Let P = U D with U and D given by 


U s 


■^ab 


1 

0 


ac 


^bc 


ha 

f^c 


"a ^ "c 


D = Diag (D-, d, D,) (3.123) 


k 
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Let e s Diag (I, m, I) and Q s dlag (0, q, 0) be dimensioned compatibly 
with (3«123) and assume P is given by the following expression. 



P = ♦ P + Q 

(3.124) 

Then the U-D 
(3.123), are 

factors of P, assumed to be partitioned consistently with 
determined as follows. 


t"ao. »o- Bol . C"ac. “ol 

(3.125) 


0^ o 

d s m d + q 

(3.126) 


“be = “ Ubc 

(3.127) 


«ab = “ ^ «ab 

(3.128) 

The matrices 

Ug and Dg satisfy the following relation. 



»a ®a "a’' = “a ■>, "a^ * q)».b "ab^ 

(3.129) 


Equations (3. 125)-(3* 128) are obtained by direct substitution of the 
partitioned U-D factors into (3.124) and by equating this expression with 
U D U^. Equation (3.129) is then easily derived with the aid of (3.126) 
and (3.128). The factorization required in Eq. (3.129) may be obtained by 
applying the Agee-Turner triangular factorization algorithm (Appendix C). 
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This single component update algorithm is the cornerstone of 
the U-D colored noise updating procedure. 


Let P = U D 0*^ where U and D are determined by equations (3.120)- 
(3.122) for the deterministic time update, (3.115). The U-D covariance 
factors after the final propagation (3.116) can be computed as follows. 


Let n and k denote the dimensions of x and p respectively. For 
I = 1,2,...,k, compute recursively equations (B.130)-(3.133) . 


*^n+l. " ^n+l 


(3.130) 


^i = '*i,n+t 


'^i.n+t = “l 


'^n+t ( 

f" ) 

^n+t / 


I i = 1,2,... ,n+l-1 


(3.131) 


'^n+l,J = H '^n+t.j 


j = n+t+1 , . . . ,n + k 


(3.132) 


Use the Agee-Turner triangular factorization (Appendix C) to compute 
the U-D factors of (3.133). 


„(A) gU) y(t)T 0 ( 1 ) g(A) 0 (i)T + c. V vT ( 3 . 133 ) 


6i 
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I 

i 

I 

t 

I 

1 

I 

J 

I 


where and 5(1) 

denote the upper left (n+t-1) rows and columns 
of the U, D matrices and 



This result is an immediate consequence of the single component 
time update, equations (3.125)-(3.129). Notice that (3.130)-(3.133) 
correspond to (3. 126)-(3.129) respectively. 


This one-at-a-time procedure also applies to the general class 
of problems involving bias parameters. If the state vector is parti- 
tioned as 


X 

P 

y 


with bias parameters in the lowest portion, then the algorithms of this 
section may be applied mutatls-mutandis. If y has dimension b then 
Eq. (3.132) must be dimensioned so that all b columns corresponding 
to y are properly scaled; i.e., Eq. (3.132) would be replaced by 


Vm = ®l^n+)l,J J = n+U1,...,n+k+b 

3.6.2 Square Root Covariance .Colored N oise Time Update 

Suppose the state error covariance for the system defined by 
Eq. (3.112) is factored as P = SS^. If S is upper triangular it may 


i 


A 
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be propagated by applying a simple modification of the mapping procedure 
developed in the previous section. The triangular square root colored 
noise update is sumnarlzed in the following algorithm. 


Triangular Square Root Colored Noise Updating Algorithm 


T 

Let P = SS* vrtjere S is partitioned such that 
n k 



(3.134) 


The propagation of S via Eq. (3.112) may be performed in two phases 
as follows. A deterministic map, Eq. (3.115), yields the intermediate 
covariance factor S, assumed to be partitioned consistently with (3.134), 
where 


Sp = Sp (3.135) 

=xp = *x =xp ♦ »xp Sp <3. 136) 

Vx’^ - («,Sj) (♦xSj)’’ (3.137) 
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The upper triangular factor Sj^ may be computed by applying the 
Schmidt time update algorithm Eqs. (3. 103)-(3. 1 1 1 ) » to the array 
W = Sx- 


The final map, Bq. (3>116), yields the covariance square root, 


S, as follows. 


For I = 1,2,...,k compute recursively Eqs. (3. 138)-(3. 144) . 


j s n+i 


(3.138) 


3.. = /tm,s,J^ + q. 


(3.139) 


Vi = S^j 




i = 1,2,...,j-1 


(3.140) 


(3.141) 


®Ji * “l®Ji 


i s J+1 , . . . ,n+k 


(3.142) 


°l= 


(3.143) 


The MGS algorithm. Section 3.3, or the Givens square .'oot factorization, 
Section 3.4, could also be used to compute Sjj. However, the Schmidt 
formula is preferred because it is more efficient t.an the other methods 
(of Chapter IV). 




rT»,=sfk%'f' .-T’ 
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Apply the Agee-Turner square root factorization (Appendix C) to 

Eq. (3.144) where denotes the upper left (n+l-1) rows and columns 

of the array S. 


g(l) g( A)T g( a) g(t)T ^ Q yyT (3*144) 

A 

This algorithm follows directly from the U-D colored noise time 
update formula by applying the identity S = 

T 

= SS where S is the general matrix 

n (3.145) 

k 

Covariance factorizations of this type are computed, for example, by 
the Potter measurement update algorithm, Eqs. (2.33)-(2.38) . Note 
that the presence of S,.„ ^ 0 introduces considerably more computation 

px 

into the propagation of S via Eq. (3.112). In order to apply the one- 
at-a^time mapping procedure, Eqs. (3. 138)- (3. 144) , it is first necessary 
to triangularize the full (n+k) x (n+k) array 


Consider the factorization P 


S = 


S S 

X ^xp 

Spx Sp 


S' = 


*x®x‘*‘*xp^px I *x^xp'*'*xp®p 


'px 


} n 
} k 


(3.146) 
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corresponding to the deterministic map in Eq. (3.115).^ For this reason, 
the one-at-a-time procedure is not an efficient method for propagating 
general covariance factors. 


Let the two-dimensional state transition be defined by 


* = I B 


and Q = 1 


If the a posteriori covariance p is propagated according to Sqs. (3.12) 
and (3.147) then 


F = P + 


Suppose P is factored such that 


A A Afi» 

p = s 


or P = U D 0' 


where 


or U s 


s diag (1,a2) 


Assume o » 1 and c? + 1 s 


'^his observation implies that propagation of the Potter square root 
for noisy systems with large numbers of bias parameters is considerably 
more costly than the corresponding U-D and triangular square root bias- 
partitioned updates (cf Eqs. (3. 120)-(3. 122) and ( 3. l35)-(3. 137)) . 
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Table 3«1 shows how the various propagation algorithms compute the 
a priori covariance factor S (or U-D). Notice that all methods except 
the RSS formula yield an S (or U-D) equivalent to t^’e rounded exact 
answer. The RSS fort.ula, however, suffers a loss of accuracy 
by applying Eq. (3.148) directly and computes a singular S with Si-| = 0. 
Slmllarlv the U-D update yields d^ s 0. 

This simple example Illustrates how use of the RSS propagation 
scheme can result In significant loss of filtering accuracy and possible 
failure of the algorithm. The other propagation algorithms, however, are 
known to be numerically reliable (cf sections 3 >3-3. 6) and thus are 
preferred over the RSS method. 
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Table 3*1< Comparison of Solutions to Examole Propagation Problem 


Exact Answer 
(Rounded) 


“ /2 o' 

1 

1 

II 

tty 

. 0 0 . 

1 — 
o 


5 = ( 2 , 0 ^) 


RSS 


U o 

0 a 




, D = (0,a^) 




1 

■ 1 r 


Givens 

- 

S = 

1 ^ = 


, D 


L 0 

c J 

Lo ij 


Schmidt 


-a 1 



(Householder) 

- 

s = 




L 0 

-0 J 




( 'i 

i 

i 


1 



■ /I 

1 

r 1 1 ^ 

f 

& 

MGS 




U = , 5 = (2,a^) 


f 



, 0 

0 J 

Lo ij 



I 5 
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Chapter IV. Cost Comparisons of Filter Algorithms 
1 Basis for Comparison 

In this chapter the various algorithms described in Chapters 
II and III are compared in terms of computer costs. Comparisons are 
based upon the arithmetic operations required by each method. Computer- 
related costs, such as indexing, storage transfers and input-output 
operations, are difficult t. j-mttfy and are neglected. However, 
these costs are c^mr-'n to a", the algorithms being considered. 

Various arithmetic operations require different execution times, 
and these differences are computer dependent. The cost comparisons 
considered here are based primarily upon the approximate execution 
times of th^ t Pacxard HP21MX-20 computer. Ti. i.s computer is 

to be . '.3 the NAVSTAR Global Positioning System '[GPS) where it 
will be employed for on-board aircraft navigation (cf General Dynamics- 
Electronics [1975]). The HP couiputer, using single precision.^ floating 
point arithmetic, has the following approximate operation times. 

■ Tj. = 3^ ijsec Tjj = 57 pSec 

Ti = 60 psec T^ = 2H00 ysec 


^Single precision arithmetic on the HP computer involves a 32 bit char- 
acteristio, or approximately 9 decimal digits. 
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Cost comparisons are slnplified if these times are normalized by the 
add time to obtain = 1.7 i Ti = 1.8 t+ , and = 70 These 
factors are used to weight the operation counts of each algorithm being 
compared.^ Thus, comparisons are based upon the relative execution 
times of the various algorithms. 

Other computers may assign somewhat different weights to the 
arithmetic operations. In most cases, however, the algorithms will 
have the same efficiency ratings, although cost differentials may vary 
from one computer to another. This point is illustrated by Including 
additional cost comparisons based upon the approximate execution times 
of the UNIVAC 1108 and Litton 4516 computers. The UNIVAC model is 
used in large, ground-based systems; while the Litton 4516 is a small, 
fixed-point computer, typical of those found on board many commercial 
and military aircraft. Hence, this cost analysis corresponds to a 
variety of computing situations. The weights associated with each 
computer involved in the study are compared in Table 4.1. 

Optimum efficiency and accuracy were the main criteria for selecting 
the computer implementation of each algorithm to be compared. Itemized 
costs for most of the mechanizations are available in Blerman [1973] 
and [1976b] and are not repeated here. However, detailed operation 
counts for the U-D algorithms may be found in Appendix E, along with 
a discussion of the colored noise propagation schemes. 


^This weighting scheme was originally suggested by Blerman [1976a]. 
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Table 4,1. Comparison of Computer Arithmetic Operation Times 


Computer 

T+ (iisec) 


Ti/ T+ 
• • 

T/-/T + 

HP21MX-20 
(32 bit) 

34 

1.7 

1.8 

70 

Litton 4516 
(32 bit) 

4 

2,2 

3.4 

250 

ONI VAC 1108 
(27 bit) 

1.9 

1.4 

4.5 

21.4 


Measurement and time updating costs are considered separately in 
order to carefully evaluate the efficiency of each individual algorithm. 
Several filter algorithms are then compared in terms of total costs at 
each stage of the filtering process. 

4.2 Measurement Update Cost Comparisons 

ine arithmetic operations required by each of the measurement 
update algorithms studied in Chapter II are indicated in Table 4.2. 

In this table n uenotes the filter dimension, while m represents the 
nuiiber of measurements included prior to variance calculations. 

The relative execution times for the HP compute.’ (cf Table 4.1) 
were used to weight the operations in Table 4.2. These weighted counts 
are included in Table 4.3 and indicate the relative execution tiroes 
of each algorithm as a function of n and m. 
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Table 4.2. Operation Counts to Process m Scalar Measurements 


Algorithm 

Adds 

Multiplies 

Divides 

Square Roots 

Conventional Kalman 






{1. + 3, 5n)m 

(!♦ 5n^ + 4. 5n)m 

m 

U 

P = P-Kl’^P 






(1 . 5n^ + 1 , 5n)m 

(1. 5n^ + 5. 5n)m 



U-D Covariance 
Factorization 

+ 

+ 

nm 

0 


(0. 5n^ - 0, 5n)* 

(n^ - n)* 



Carlson Triangular 
Covariance 
Square Root 

(1. 5n^ + 2. 5n)m 
+ 

(0. 5n^ + 0. 5n)* 

(2n^ + 5n)m 
+ 

(0.5n^ + 0. 5n)*^ 

(2n 4* l)m 

nm 


(3n^ H 3n)rn 

{3n^ + 4n)m 



Potter Covariance 
Square Root 

f 

+ 

2m 

m 






Stabilized Kalman 





P = (I-Ka'^) P(l-KaT)'r 

(4n^ 4 4n)m 

(4n^ + 6n)m 

m 

0 

+ Ktk'^ 






*These operations involve computing estimate error variances. 
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Measurement update costs are further compared in Figure 4.1. 
Execution times have been normalized by the corresponding conventional 
Kalman times and are display ad as a function of n for m = 1 and m = 5. 
Hence, this figure gives the various algorithm costs as a percentage 
increase over the conventional updating costs. Notice that the U-D 
algorithm is considerably more efficient than any other alternative 
to the conventional Kalman method. Moreover, when variances are required 
infrequently, the U-D and conventional Kalman costs are nearly identical. 


Table 4.3. Measurement Update Costs Weighted for the HP Computer 


Algorithm 

Execution Tiine/T+ 

Conventional Kalman 

(4. 1n^ + 1 1 .2n)m 

U-D 

(4.1n^ + 1?.7n)m + (2.2n^ - 2.2n)» 

Carlson 

(4.9n^ + 84. 6n + 1.8)m + (1.4n^ + 1.4n)* 

Potter 

(8.1n^ + Q.8n + 73.6)m + (2.7n^)* 

Stabilized Kalman 

(10. 8n^ + I4.2n)m 

•Variances computed 


In sharp contrast to this efficiency, the other .'.Igorithrns generally 
require at least 100!^ more computation time than the conventional method. 
Note that the Carlson normalized costs Increase as the state dimension 
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Figure 4.1. Measurement Update Cost Comparison 


decreases. This trend is related to the n scalar square root calculations 
in the Carlson formula; these operations are relatively more expensive 
when n is small. 

4.3 Time Update Cost Comparisons 

4.3.1 Comparison of General Propagation Algorithms 

Several algorithms for propagating covariance factors have 
been described in Chapter III. Wi-h the exception of the RSS methods, 
all of the algorithms are numerically reliable. Therefore, the most 
appropriate time update algorithm for each covariance factorization 
may be selected on the basis of efficiency. 
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Consider the covariance square root propagation schemes 
listed in Table 4.4. The operation counts given for each algorithm 
include only the calculations required to compute the updated square 
root S from the n x (n+k) array W = i ♦§]. Notice from Table 

4.4 that the Schmidt and RSS methods require less computation than 
the other methods. Note, also, that these two algorithms involve equiv- 
alent amounts of computation. Therefore, no loss of efficiency is 
encountered when the more reliable Schmidt algorithm is selected for 
covariance square root propagation. The Schmidt time update may be 
coupled with the Carlson or the Potter square root data processing 
algorithms. The resulting filters are referred to as the Carlson-Schmidt 


Table 4.4. Arithmetic Operations Required by Square Root 
Triangularization Algorithms 



^ i 

' \ 

' i 75 


I 




and Potter-Schmidt algorithms. These two square roc-v covariance filters 
will be considered in the cost comparisons to follow. 

A similar comparison of U-D propagation methods is given 
in Table 4.5. The calculations listen for each algorithm involve only 
the computation of U and D from the arrays W = [B I ♦!)] and 
D = diag (Q, D). The modified Givens operation counts in Table 4.5 
represent the minimum attainable with this method; l.e., the more effi- 
cient formula, Method B in Eq. (3.87), is applied a maximum number of 
times. This assumption is justified by the final dijcuscion in Section 
3.4. The reader may wish to refer to Appendix E for a more detailed 
description of the U-D propagation costs. 

It is not immediately apparent from Table 4.5 which U-D 
algorithm is the most efficient. However, it is clear that the RSS 
method has very little cost advantage over the modified Givens algorithm. 
Of the two reliable U-D methods, the Givens method is more efficient 
when n is large. However, when n is small and k is large the extra 
divide operations required by the Givens method are more apparent, 
and consequently the MWGS algorithm may be the more efficient method. 
Hence, in the cost comparisons to follow both of these U-D propai^jation 
schemes are considered. 


REPRODUCIBILrrY OF THE 
ORIGINAL PAGE IS POOR 
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Table 4.5. Arithmetic Operations Required By U-D 
Factorization Algorithms 


Algorithm 

Adds 

Multiplies 

Divides 

Square Roots 


0,7n^ - 0.5n^ + 0.8n 

0.7n^ + 2.5n^ - 0.2n 



Modified Givene 
(Eqe. (3.80)-(3.89)) 

+ 

+ 

0. 5n^ - 0, 5n + nk 

— 


n k 

(n* + 4n)k 





+ 1.5r^ . 0.5n 



MWGS 

(Eqt. (3,28)-(3e30)) 

3 ’ 

n + n^k 

+ 

n - 1 

— 



(n** + n)k 




0. 7n^ + + 0. 3n 

0. 7n^ + 2. 5n^ -l.2n 


i 

RSS 

(Section 3,2) 

+ 

+ 

n - 1 

1 _ 

1 


(0,5n^ + 0,5n)k 

(0.5n^ + i.5n)k 




The arithmetic operations involved in a general time update 
via Eq. (3.1) are listed in Table 4.6 for each of the filter algorithms 
to be considered. These counts include the calculations required to 
compute estimates and the necessary updated covariance factors. The 
differences between the Potter-Schmidt and Carlson-Schmidt propagation 
costs are due entirely to the different structures of the covariance 
square roots involved. For example, computation of the product ♦S 
with a Potter square root requires 0.5n^ more adds and multiplies than 
the same calculation involving a Carlson triangular square root. 

Table 4.7 contains the weighted propagation costs for 
each filter. These costs represent the relative execution times of 
the various time update schemes when the HP computer is used. Costs 
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are further compared in Figures 4.2 and 4.3 as a function of k/n where 
k denotes the number of process noise parameters. Execution times 
have been normalized by the corresponding execution time for the con- 
ventional method. Hence these figures show the percentage of conventional 
mapping costs which are required by the different algorithms. The costs 
associated with variance computations in the factorization filters 
are omitted from Figures 4.2, 4.3 and those to follow. However, the 
relative algorithm costs are not altered appreciably by Including these 
computations (cf Table 4.7). In fact, when variance calculations are 
Included, the cost curves in Figures 4.2 and 4.3 experience an upward 
shift of at most three percent. 

Notice in Figures 4.2 and 4.3 that normalized propagation 
coSwS Increase as a function of k/n. However, even when k/n = 2 all of 
the algorithms require less than 50 % more computation than the conven- 
tional propagation formula. The Potter-Schmidt update is the most 
expensive time update algorithm while the other methods are generally 
competitive with one another. The modified Givens method is typically 
the more efficient U-D propagation scheme although for small dimensioned 
systems with many process noise parameters (l.e., k/n > 1), the MWGS 
algorithm may be less expensive (cf Figure 4.2). In either figure the 
cost differentials between the MWGS and the Givens methods are modest, 
and both schemes are generally competitive with the conventional propa- 
gation algorithm. Subsequent discussions of U-D propagation methods 
are therefore limited to the modified Givens algorithm. 
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for a General Time Update 


Algorithm 


Conventional 

Kalman 

1. 5n^ 4 2n^ 4 0. 5n 
4 

(O.en^ 4 0. 5n)k 


l,5n^ 4 0.5n^ 


4 

U-D (MWCS) 

n^k 


4 


(0.5n^ - 0.5n)* 


U-D (Givens) 


1.2n' 4 0.8n 
+ 


Multiplies 


1.5n^ + K5n^ 
+ 

(0.5n^ + 1.5n)k 


1.5n^ 4 2n^ - 0. 5n 


(n" -, ^)k 


Divides I Square Roots 


(n" - n)’> 


1.2n^ 4 - 0. 2n 


in 4 4n)k 
4 


(n^ - n)« 


- 0. 5n 
4 


Carlson- 

Schmidt 


1. 2n 4 1 , 5n h 1, 3n 
4 


1. 2n 4 Z , 5n 4 0. 3n 
4 

(n^ 4 n)k 
4 


(0., a" 4 0,5n)^ 


(0, 5n^ 4 0. 5n)* 


1.7n^ 4 4 1.3n 


1.7n^ 4 2n^ 4 0,3n 


Potter- 

Schmidt 


(n'^ 4 n)k 


♦Variances computed. 
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Table 4.7. General Time Update Costs ! .sighted for the HP Computer 


Algorithm 

Execution Tlme/t^. 

Conventional Kalman 

4.1n^ + 4.6n^ + 0.5n + (1.4n^ ♦ 3.1n)k 


4.1n^ + 4.8n^ + O.ln + (2.7n^ + 1.7n)k 

U-D fMWGS) 

+ 


(2.2n^ - 2.2n)» 


3.2n^ + 6n^ - .4n + (2.7n^ ♦ 8.6n)k 

U-D (Givens) 

+ 


(2.2n^ - 2.2n)« 


3.2n3 + 5.8n^ + 73. 6n + (2.7n^ + 1.7n)k 

Car] son-Schmidt 

+ 


(1.4n2 + 1.4n)» 


4.6n3 + 4.4n^ + 73-6n + ( 2 . 70 ^ + 1.7n)k 
Potter-Schmidt + 


(2.7n^ + 2.7n)« 


*Varlances computed 
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k/n 

Figure 4.2. General Tine Update Cost Comparison, n = 10 



Figure 4.3. Gener<^l Time Update Cost Comparison, ti s 30 
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4.3.2 Comparison of Colored Noise Propagation Algorithms 

The operation counts associated with a colored noise time 

A 

upuate of the various covariancs factors are compared in Table 4.8. 

The U-D algorithm to be considered propagates process noise one component 
at a time and employs a modified Givens factorization for the determin- 
istic phase of the update (of Eqs. (3. 120)-(3. 122) and (3.130)-(3. 133)). 

A similar algorithm based upo. ;i.e MWGS factorization requires approxi- 
mately the same computation time as the one considered here (see Appen- 
dix E). 


Colored noise time updating of the Potter and Carlson 
covariance square r''.'t:; is mor*^ efficiently accomplished by mapping 
the process noise in one step. Hence, the Potter-Sehmidt and Carlson- 
Schmidt operation counts in Table 4.8 represent the cost of a Schmidt 
time update which exploits system structure. The reader is referred 
to Appendix E for further details of these mechanizations. 

The weighted operation counts in Table 4.9 repre'».ent the 
relative costs of a colored noise update on the HP computer. These 
costs, normalized by the conventional propagation costs, are iUuc‘”'eied 
in Figures 4.4 and 4.5 as a function of k/n. Filter dimension in these 
figures is n+k where k represents the numb- of colored noise parameters. 
For all values of k and n the ®otter-Schmidt. colored ncise update is 
considerably less efficient thur the Carlson-Schmidt and U-D algorlchms. 

^The colored noise system is deflnsd by E. 7 . ' 3. (12). 
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Table 4.8. Operation Counts for Colored Noise Time Update 


) Algorithr 


Conventional 

Kalman 


Carlson- 

Schmidt 


Potter - 
Schmidt 


Adds 

Multiplies 

1. 5n^ + 1. St? + nk^ 

1, Sn^ + i. 5n^ + (n + l)k^ 

+ (2. 5n^ +1 Sn e l)k 

+ (2. 5n^ + 2. 5n + : /k 

1.2n^ + 0.8n + 0, 3k^ 

1. 2n^ + 3n^ - 0, 2n 

+ (1..5n - 0.5)k^ 

+ 0,3k^ + (1.5n + 2)k^ 

+ (Zt? - 0. 5n + 1. 2)k 

+ (2n^ + 3n + 0. 7)k 

+ (0. 5(n + k)^)* 

+ {(n + k)^ - (n + k))* 

- 0. 5(n + k))* 


1.2n^ + I.5n^ + 1.3n 

1. 2n^ + 2. 5n^ + 0. 3n 

+ 0. 3k^ 

+ 0. 3k^ 

+ (l,5n + l.'jk^ 

+ (1.5n + 2.5)k^ 

+ (2n^ + 4. Sn + 1. 2)k 

+ (2n^ + 5. 5n + 2. 2)k 

+ (0. 5(n + k)^)* 

+ (0. 5(n + k)^)* 

+ 0. 5(n + k))* 

+ 0. 5(n + k))* 

1.7n^ + + 1.3n + k^ 

1.7n^ T 2n^ + 0. 3n + k^ 

+ (n + 2. 5)k^ 

+ (n + 4)k^ 

+ (3n* + 3n + 2. 5)k 

+ (3n^ + 4n + 3)k 

2* 

+ (n + k)^ 

2* 

+ (n + kr 


Divide. 


n + k - 1 


* Variances computed 


























33-798 


Table 

4.9. Colored Noise Time Update Costs 
Weighted for the HP Computer 

Algorithm 

Execution Tlme/x+ 

Conventional 

Kalman 

4.1n3 + 4.1n^ + (2.7n + 1.7)k^ + (6.8n^ + 5.8n + 4.4)k 


3.2n3 ♦ 6n^ - .4n + 0.8k^ + (4. In + 3.8)k^ 

U-D 

^ + 


(5.4n^ + 6.4n + 3.3)k + (2.2(n + k)^ - 2.2(n + k))» 


3.2n3 + 5.8n2 + 73. 7n + 0.8p3 + ( 4 . in + 5.8)k^ 

Carlson- 

Schmldt 

+ 

(5.4n2 + 13 . 9n ♦ 76.7)k + (1.4(n + k)^ + 1.4(n + k))» 


32 , 5 

4.6n + 4.4n + 73. 7n + 2.7k^ + (2.7n + 9.3)k‘^ 

Potter- 

Schmidt 

+ 

(8.1n^ + 8. In + 79.4)k + (2.7(n + k)^ + 2.7(n + k))» 

*Varlances computed 


Unlike these triangular factorization methods, the Potter square root 
does not allow for full exploitation of the special system structure 
In Eq. (3.112). The Carlson-Schmldt and U-D colored noise algorithms 
are generally competitive with one another, although the U-D method 
Is somewhat more efficient when ns 10. 


Note that cost differentials of all the methods Increase 


as a function of k/n. owever, even when k/n s 2 the U-D costs are 
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within 30jt of those related to the conventional propagation method. 

Moreover, when k/n < 1 the U-D method is the least expensive colored 
noise algorithm. 

Of the two figures being discussed, Figure 4.4 is more repre- 
sentative of the costs usually incurred with small, on-board computers 
of the type considered here. For example, the GPS project expects to 
employ the HP computer for airborne filtering of a colored noise system 
which has dimensions n = k s 6 (cf General Dynamics-Electronics [1975]). 

Higher dimension filters, such as those represented by Figure 4.5, 
are generally designed for ground-based analysis on large computers 
such as the UNIVAC 1108 or IBM 360 s Thus, further comparisons of colored 
noise costs will be restricted to the case n = 10. 

£ = 

4.4 Comparison of Total Filtering Costs 

The measurement and time update costs given in the two previous 
sections may be combined to yield a rough comparison of total filtering 
costs. Most of these comparisons are based upon the assumption that 
a single scalar measurement is Included at each step of the filtering 
process. Note that measurement updating is generally an order of magni- 
tude less expensive than time updating (cf Tables 4.3 and 4.7). Hence, 
filtering costs primarily reflect time update expenses, especially when 
n is large. For small systems, however, measurement update costs are 
more visible, and the addition of several observations may noticeably 
alter the filter cost differentials. To Illustrate these effects- . 
additional comparisons involving multiple measurements are included 
for the case n s 10. 
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Figure 4.4. Colored Noise Time Update Cost 
Comparison, n = 10 



Figure 4.5. Colored Noise Time Update Cost 
Comparison, n s 30 


reproducibiuty op ran 
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Filter cost comparisons do not Include the expenses related to 
variance calculations. As previously noted, these computations have 
a modest effect on relative time update costs^ and, hence, on tue normal- 
ized filter costs to be compared. 

Filtering costs associated with the general system model. Eg. (3*1), 
are compared in Figures 4.6 and 4.7. These figures represent the relative 
computation times required by each filter algorlthmi^'*^ and were obtained 
by combining the weighted execution times In Tables 4.3 and 4.7. The 
costs of each algorithm have been normalized by the corresponding conven- 
tional Kalman costs and are displayed as a function of k/n for n s 10 
and n = 30. In a similar manner, the colored noise filter algorithms 
are compared in Figure 4.8. Each colored noise filter has dimension 
n+k where n = 10. 

Notice first that cost differentials associated with the covariance 
factorization methods Increase as a function of k/n, while those related 
to the stabilized Kalman algorithm remain relatively constant. In 
each case the stabilized Kalman costs are within 20f of those related 
to the conventional Kalman formula. Of all the algorithms compared, 
the Potter-Sohmldt method is the most expensive. The Potter-Schmldt 
colored noise algorithm is particularly Inefficient when there are 


^See the discussion of Figures 4.2 and 4.3. 

^^The U-D filter algorithm included in these comparisons employs the 
modified Oivens time update method. Costs associated with this filter 
algorithm may be easily adjusted to obtain comparisons for a U-D filter 
based upon the MWOS time update algorli-hm (of Figures 4.2 and 4.3). 
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large numbers of colored noise parameters. For example, when k/n s 2 
this algorithm is tw'.oe as expensive as the conventional Kalman method. 

By oomparison, the Carlson-Sohmidt and U-D colored noise algorithms 
are significantly more efficient and require less than 50 % additional 
computation relative to the conventional Kalman formula. 

Notice that the U-D method is the least expensive covariance 
factorization algorithm. When system dimension is small, this method 
has a noticeable cost advantage (cf Figures 4.6 and 4.8). Moreover, 

U-D costs are generally within 20 % of those related to the conventional 
Kalman method, and when k/n < 0.4 the U-D algorithm is the least expensive 
method considered . 

For large scale systems there Is less diversity in filter algorithm 
costs than that observed in Figures 4.6 and 4.8 where ns 10. Consider 
the comparison in Figure 4.7 for the case n = 30. Notice that filter 
algorithm costs differ less in this comparison than they do in the corre- 
sponding lower dimension case in Figure 4.6. Cost differences are less 
apparent in Figure 4.7 because time updating expenses dominate the compari- 
sons when n s 30. Compare Figures 4.3 and 4.7 and note that relative 
filtering costs are nearly identical to relative propagation costs. 

Although filtering costs primarily reflect the expenses related to 
time updating, the cost effects of a single measurement update are apparent 
when n s 10. Compare, for example, Figures 4.2 and 4.6. Notice that 
the U-D results in each figure are nearly identical, but the Potter- 
Sohmidt and Carlson-Schmidt cost curves are noticeably different. These 
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PERCENTAGE 
INCREASE 
: OVER 

CONVENTIONAL 
KALMAN 
COSTS 



0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 


k/n 


Figure 4.6. Cost Comparison of General Filter 
Algorithms, n s 10 


PERCENTAGE 
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OVER 
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KALMAN 
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k/n 


Figure 4.7. Cost Comparison of General Filter 
Algorithms, n a 30 
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Figure 4.8. Cost Comparison of Colored Noise Filter 
Algorithms, n s 10 
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differences reflect the measurement cost differentials in Figure 4.1 
where only the 0-D algorithm has costs comparable to the conventional 
Kalman isethod. Thus, we would expect additional measurement updates 
to Increase all of the normalized filter costs in Figure 4.6 except 
those corresponding to the U-D algorithm. 

Consider the case where three measurements are Incorporated 
after each time update. Relative filtering costs for this case are 
compared in Figures 4.9 and 4.10. The results in Figure ^(.9 correspond 
to the general system model, while Figure. 4.10 illustrates the colored 
noise filter costs. By comparing these figures with the corresponding 
results in Figures 4.6 and 4.8 we find that, except for the U-D method, 
each filter algorithm becomes relatively more expensive as additional 
measurements are incorporated. Figures 4.9 and 4.10 indicate that 
the Carlson-Schmidt and Potter-Schmidt lower-dimensioned filters can 
have equivalent costs when multiple measurements are involved. Further 
note from Figure 4.10 that whenever k/n < 1 the U-D method is the least 
expensive colored noise filter algorithm and requires at least 4oi( 
less computation than the other covariance factorization methods. 

The cost comparisons discussed in this chapter have all been related 
to the HP computer. The relative costs of the various algorithms may 
vary considerably, depending upon which computer is used. Consider 
the different weights that the Litton and UNIVAC cenputers asalgn to 
the various arithmetic operations (see Table 4.1). Note that square 
root calculations are ten times more expensive on the Litton computer 
than on the UNIVAC model. Hence, we would expect the Carlson-Sohmidt 


91 







33-798 

and Potter-Schmidt covariance square-root methods to be noticeably 
less efficient on the Litton computer. 


Cost comparisons related to the Litton and UPIVAC computers are 
Illustrated In Figures 4.11-4.14 for the case where n s 10. These 
oom^rlsons were obtained In the same manner used to generate the costs 
In Figures 4 and 4^ weighting facers differ. 



Notice from Figures 4.12 and 4.14 that when the UNIVAC computer 
Is used, algorithm cost dlfiferenoes are generally modest. However, 
filtering costs related to the Litton model differ considerably. On 
the Litton computer the Carlson-Sohmldt and Potter-Schmidt colored 
noise algorithms have nearly equivalent costs and are generally 60-80K 
more expensive than the corresponding U-D method. By comparing Figures 
4.11-4.14 with the corresponding HP costs in Figures 4.6 and 4.8 we 
find that the relative costs of the U-D method are approximately the 
same on each computer involved In this analysis. In each case 0-D costs 
are typically within 20$ of those related to the conventional Kalman 
method, and when k/n <0.4 the U-D algorithm Is the most efficient 



method studied . 

4.5 gonolufllona 

Cost analysis has demonstrated the relative efficiency of the U-D 
filtering technique. This method, employing either the modified Givens 
or MHOS propagation algorithms. Is generally competitive with the oonven- 
tional Kalman formula. The U-D filter algorithm based upon the modified 
Givens time update la particularly efficient and for problems Involving 
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Figure 4.9. Relative Costs of General Filter Algorithms 
When Three Measurements are Included, n s 10 
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Figure 4.10, Relative Costa of Colored Noise Filter 
Algorithms When Three Measurements are Included, 
n « 10 


» gi 








Figure 4.11. Pelative Coste of General Filter Algorithms 
When Litton Computer is Used, n = 10 



Figure 4.12. Relative Costa of General Filter Algorithms 
tlhen UNIVAC Computer is Used, n s 10 
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Figure A. 13 . Relative Costs of Colored Noise Filter 
Algorithms When Litton Computer is Used, n s 10 
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Figure A. 14. Relative Costs of Colored Noise Filter 
Algorithms When UNIVAC Computer is Used, n s 10 
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modest numbers of process noise parameters, this U-D algorithm is the 
least expensive of all the methods considered. By comparison, the 
Potter-Schmldt and Carlson-Schmidt square root algorithms are signifi- 
cantly more costly, particularly in real-time appliOations involving 
small computers. In these situations the Carlson-Schmidt method is 
generally no more efficient than the Potter-Schmldt algorithm and can 
be more expensive. Moreover, both methods can involve in excess of 
60 % more c(xnputation than the corresponding U-D algorithm. 
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Chapter V. Covariance Error Analysis Techniques 


5.1 Introduction 

Several difficulties are frequently encountered when filtering 
algorithms are applied to an actual estimation problem. For example, 
precise knowledge of the system model and complete a priori statistics 
are often unavailable. In addition, computational constraints may 
limit the dimensionality of the filter. Because of these problems 
suboptlmal filtering Is often Inevitable. As an Important aid In filter 
design, one may perform a sensitivity analysis of the estimate error 
covariance. Techniques for evaluating filtering sensitivity to Incorrect 
a priori statistics have been derived for the discrete case by Fagln 
[1964], Nlshlmura [1966 and 1967] and Heffes [1966]. Their results 
were extended by Griffin and Sage [ 1968] to Include the effects of 
Incorrect data and state transition matrices. These error analysis 
methods all propagate an actual, or true, error covariance. The algo- 
rithms are often cumbersome and computationally expensive. Furthermore 
their reliability Is questionable since they employ the numerically 
unstable Kalman formulas. 


t 


In this chapter a new approach to covariance error analysis Is :■ 

presented which Is based upon the U-D covariance factorization. This j 

method facilitates an accurate evaluation of general modeling errors, 
including a) incorrect a priori statistics, b) mismodeled system dynamics 
and data equations, and c) Incomplete parameter sets. A general error 
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analysis algorithm which Incorporates all of these effects is derived 
in two stages. 

The case of incorrect a priori statistics is considered first. 

This case includes a large class of modeling errors which are easily 

-«■ 

analyzed. The basic evaluation algorithm developed for this Important 
subproblem is then extended to include the effects of Incorrect data 
and state transition matrices. The extended analysis provides a general 
evaluation algorithm capable of considering the effects of unmodeled 
or neglected parameters. However, the analysis of unmodeled bias param- 
eters is included as a separate topic in order to highlight the simplicity 
and flexibility of the resulting algorithm. 

5.2 Evaluation of Incorrect A Priori Statistics 

Suppose a discrete linear filter is designed according to the 
following model 


x(i+1) = *(i)x(i) + B(i)w(i)) (5.1) 

> i=0 12 * * * 

z(i) = A(i)x(i) + v(i) ) ’ ' ’ (5.2)^ 


where 


X e ; w e Rj^ ; z e R^ 


t 

Scalar measurement coefficient matrices are denoted as A Instead of 
the "a^" used in Chapter II . For the error analysis discussed in this 
chapter the symbol; "is more appropriate, particularly in section 5.2 
where the notation A^ is required. 
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x(0) - N(0,P) ; w(i) - N(0,Q(i)); v(i) - N(0,r(i)) (5.3) 

and 

E{x(0)w(i)^)=0; E{x(0)v(i)}s0; E{w(i)v( j)}=0; i,J>0 (5.4) 


Let the correct process and data model be the following 


Xj^(i+1) = *(i)x3(i) + 
z(i) = A(i)x^(i) + VgU) 


i— 0 |1j2y... 


(5.5) 

(5.6) 


where 


x^(0) 


Xa(i), Wg(i) and v^(i) satisfy (5.4) and 
N(0,Pa(0)) ; w^(i) - N(0,Qa(i)) ; V^(i) - N(0,ra(i)) 


(5.7)^ 


Let G(i) denote the computed (suboptimal) gain at stage i. 

Then the actual estimate errors associated with the model (5.5)-(5.7) 
satisfy the following recursions 


A^(i) = x(i)-x^(i) = (I-G(i)A(i))Ax(i) + G(i)v^(i) (5.8) 

A5l(i+1) = 3f(i+1)-x^(i+1) = *(i)Ax(i)-Bg(i)Wj^(i) (5.9) 

where 

A5K0) = -Xa(0) (5.10) 

Equations (5.8)-(5.10) may be verified by applying (3-11). (3-13) and 
(3.16) together with (5.5)-(5.7). 


There is no loss of generality in assuming a zero mean process. This 
assumption is made in order to free the error analysis of superfluous 
algebra . 


99 


33-798 


Thus, the actual a priori and a posteriori error covariance for 
this problem may be computed from the following recursions 

= ( I-GA )P3( I-GA + Gr^G''^ (5.11) 

\ + Wa^ <5.12) 

where initially 

\ s P^CO) (5.13) 

Equations (5.11)-(5.13) represent the conventional error analysis 
algorithm for mlsmodeled a priori statistics (cf Nishimura [1966] 
and Heffes [1966]). This evaluation method is simply an application 
of the stabilized Kalman covariance equations and for this reason it 
is numerically unreliable.^ 

A more reliable error analysis method may be obtained by employing 
the U-D factorization techniques studied in Chapters II and III. Time 
propagation of the actual U-D covariance factors, via (5.12), may be 
achieved by a direct application of the MWGS algorithm (or the modified 
Givens update) in Chapter III. The U-D suboptimal measurement update, 
corresponding to (5.11), is based upon the following covariance decomposi- 
tion. 


tThe numerical deficiencies of the stabilized Kalman formula are demon- 
strated in Chapter VI where, in sharp contrast, the U-D and Potter 
methods are shown to possess superior numerical characteristics. 
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Given: 


P - a priori error covariance 


A.r^ - observation coefficients and scalar measurement 
error covariance 
G - arbitrary gain 


The a posteriori error covariance where 


P = (I-GA)P(I-GA)^ ♦ Gr^G^ 


may be written as follows 


(5.14) 


where 


P = P + oU^ 


P = F-Kv*^ 
1 s K-G 


V = PA**’ 


(5.15) 


(5.16) 

(5.17) 

( 5 . 18 ) 


s v/a 


a = Av + r„ 


(5.19) 


Consider the quadratic P(G) in Eq. (5.14), Note that the quantity 
K defined by (5.19) minimizes P(G) and that this minimum is P(K) = F. 
Since the difference P(G)-P is non-negative definite and symmetric, 

A 

we would expect P(G) to have the form (5.15). This decomposition is 
readily verified by substituting (K-X) for G in Eq. (5.14), regrouping 
terms and applying the identities (5.18) and (5.19). The covariance 
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P thus corresponds to the minimum quadratic, i.e., the optimum updated 
covariance associated with the Kalman gain K . ^ 

Suppose the covariance P in (5.14) is factored such that 
P s UDo'*’. Equation (5.16) then implies that the U-D factors of P may 
be computed from the Bierman optimum measurement update, Eqs. (2.69)- 
(2.79). Given the 0-D factors of P, Eq. (5.15) suggests using the Agee- 
Turner triangular matrix factorization (Appendix C). Thus, the U-D 
suboptimal measurement update may be performed in two stages; an optimal 
update to a modified problem, ^ followed by an Agee-Turner rank one 
matrix update. This method is summarized in the following algorithm. 

U-D Arbitrary Gain Update Algorithm 

Given; U,D - factors of a priori covariance P 

A,r^ - observation coefficients and 

scalar measurement error covariance 
G - arbitrary gain 

The a posteriori covariance factors 0 and 6 can be obtained as follows: 

• Compute via the U-D optimum measurement update algorithm, 
Eqs. (2.69)*(2.79) , the quantities: 


^The gain K is optimum given a priori covariance F. However, F does 
not represent the optimum (or minimum) covariance attainable unless 
P is also optimum. 
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K - Kalman gain 
a ~ innovations covariance 
U,D - the factors of F s (I-KA)F 


• Set X ~ K— G 


• Apply the Agee-Turner algorithm, Appendix C, to compute 
0 and ^ where 

_ Q50T ^ (5.20) 

This evaluation may be represented symbolically as 


[U.S.A.rg] ♦ [U,D,x,o] ♦ CO.S] 


The U-D arbitrary gain update inherits the efficiency and reli- 
ability of the two algorithms it employs. The optimal U-D update was 

p 

shown in Chapter IV to require 1.5n + 0(n) multiplications while the 

2 

Agee-Turner factorization requires n such operations (cf Appendix C). 

2 

Hence the entire update is performed with 2.5n 0(n) multiplications. 

When variances are computed, total multiplications reach 3n + 0(n) 
as compared to the 4n^ 0(n) such calculations required by the conven- 

tional evaluation formula, Eq. (5.11). A detailed cost comparison 
of gain evaluation algorithms is included in Appendix F. 


This arbitrary gain update is the cornerstone of the U-D error 
analysis technique. Note that evaluation of incorrect a priori statistics 
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requires only a minor addendum to the optimal U-D filtering algorithm. 
Thus, both the gain computations and the subsequent error analysis 
can be performed by one program. The gains may be computed during 
a first pass through the filter, operating in what it believes to be 
an optimum mode. The same filter, tuned to account for suboptimal 
gains and given the correct system model, may then perform the gain 
evaluation . 

A similar error analysis algorithm could also be developed for 
the factorization P s SS'' where S is upper triangular. The suboptimal 
measurement update for this factorization closely parallels the 0-D 
formula. First, S, the square root of P in Eq. (5.15), is computed 
via Carlson's optimal update, Eqs. (2.39)-(2.46) . Then an Agee-Turner 
square-root factorization (Appendix C) yields ^ where = ss"^ + axx"^. 
This update may be coupled with a Schmidt (Householder) propaga ion 
scheme to obtain a complete square root error analysis algorithm. 

The square root method is not pursued further, however, because it 
lacks the efficiency of the U-D evaluation algorithm (see Appendix F). 

5.3 EYalufttjl9P..<;if.Mi.9Bi9d?ls0.fttfl-and-Stots,Tr.flnatU9n 
5.3.1 ggpgpfti E,r/9r APftlygt? T99hb^qw9 

In this section the U-D error analysis method- is extended 
to account for mismodeled data and state transition matrices. Suppose 
the correct process and data equations are the following 
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( 5 . 21 ) 

( 5 . 22 ) 


x^a+l) = ♦a(i)Xj^(i) + Ba(l)Wa(i) 

z(l) s Aa(i)Xa(i) + Va(i) 
where w^t and *g(8) satisfy Eq. (5.7). 

Suppose the filter applied to this problem is based upon. the 
assumptions (5.1 )**(5.4 ) and generals the j^ain seiiuenee .; Theh 

the agtu ffi * posteriori: error equation at each stage has the form 


where 


where 


aft = ^-X^ S (I-GA)AX ♦ QAAXj^ ♦ GVg 

,5.23)^ 

AA = Ag-A 

(5.24) 

r, iSi, is related to previous errors as 

follows 

iSi s x*Xg = *Aft - AfcCjj - BgWg 

(5.25) 

A* = •a - * 

(5.26) 


Due to the errors AA and A*, covariance expressions for A^ and 
ASr contain the actual state covariance, = Elx^Xa^), and the cress- 
covariance, * E{AxxJ). Thus, covariance error analysis for this 
problem J^iwelves propagation of the 2n-dimensional covariance 



+Hote that mismodeled data and state transition matrices produce biased 
estimates unless E{x,} > 0. 
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The recursions required to propagate S* by conventional methods were 
derived by Griffin and Sage [1968]. However, their algorithm appears 
unwieldy and, Judging from the experimental results in Chapter VI, 
it is numerically unsound. The problem may be solved more accurately 
and efficiently by applying 0-D covariance factorization methods. 

For this error analysis we will find it convenient to avoid covariance 
equations altogether and study. Instead, the associated "weighted error" 
propagation problem.^ To this end, we' interpret the U-D covariance 
factorization as a particular "idiltening" of the estimate errors. 

That is, in terms of second-order statistics the errors may be described 
as 


AX = Ov^ 


(5.27) 


where the zero mean noise, has diagonal covariance D. 


Suppose the U— D factors of the a priori covariance 


0m 

.0 S 




are partitioned as follom. 



n n 


diag (I^,Djj) 


(5.28) 


^Thls approach is analogous to the error analysis technique associated 
with square-root information filtering (of Bierman [1976b]). 
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Then the •’whitening" interpretation of Bq, (5.27) implies that a priori 
estimate errors may be written as 


t! = !!...!!« !! 

X- 0 U* 


(5.29) 


where v-j and are zero mean, uncorrelated and have diagonal covariances 
Dj and respectively. Equations (5.23) and (5.29) yield the following 
a posteriori error equation 


where 


' ‘ '*a 

‘I-O*)"* % -- , 

)n 

0 Ujj - 

"J V, )n 


U~ s (I-GA)U«^ + GAAU. 


(5.30) 


(5.31) 


The noise vector » v, has zero mean and diagonal covariance 


1 n n 


D . (r^, Djf, D^) 


Given the expressions (5. 30) *(5.32)1 we dosire the representation 


(5.32) 


Ufcc h»2 


r • :-t”r .’v 


M h x 


(5.33) 


107 




r 



with 


33-798 




n n 


6 S diag (Dj^, D') 


The constraint that 


(5.34) 



be unchanged during the transforma*’ on from (5.30) to (5.33i^^ yields 
the identities 



CU', D'] s [Ujj, Djj] (5.35) 

Ufbc = % “ «5bc - < 5 . 36 ) 

U-DjUjT * (I-OA)U5D3jU5f'^(I-GA)^ + GraO'^ (5.37) 

Thus, Uj^, and d' are obtained directly. The factors Uj^ and 
may be computed by applying the U-D arbitrary gain algorithm (section 
5.2) to Eq. (5.37). 

Similarly, a U-D algorithm may be derived for the time update, 

Bq. (5.25). Vhen the expression (5.33) is applied to (5.25), the a priori 
error equation takes the form 
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r T r _ T *^a 

M K ^st Hx “ , 

y=|T;[-Ty 5 ; 


(8.38) 


% = *OSx - 


(5.39) 


The noise vector v in Eq. (5.38). has zero mean and covcrianoe 5 where 


k n n 


5 = diag (Qg, Dj, \) 


(5.40) 


Note that one of the U-D time update algorithms in Chapter III may 
be used to triangularize (5.38). In this ease the W and D arrays involved 
it', the time update are given by Eqs. (5.40) and (5.41). 


k n n 


I ”8 1 % 1 1“ 
>'n I Vxj )» 


(5.41) 


Computer storage and computation required for this update may be reduced 
by taking advantage of the nxn b?ook of zeroes in the H array. 


The general U-D error analys' i method is summarized as follows. 
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General O-D Covariance Error Analysis Algorithm 


Given: {*(!)) and {A(i)) - the incorrect state transition and 

data matrices assumed by the filter 

- suboptimal gain profile computed by the filter 

{*g(i)} - correct state transition matrices 

|Bg(i)} - actual process noise coefficients 

{Ag} - actual data coefficients 

{r_(i)} - actual data noise covariances 

{Qa(i)) - actual process noise covariances 

cl 

ana u = diag (D~, - factors of the 

true initial covariance of x and x^.^ 

The true error covariance for this model is propagated according to 
the following recursions. 



^Typically, the initial factors ^re such that 
[Ujf, Djf] = tUjj, Djj] and * 0, whore 


P,(0). 


BBPRQDOCIBn^ OB^ 
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= ^*x ” °^*%x ■ AAO^)} AA r Ag - a (5.112) 

Apply the Bierman algorithm, Eqs. (2. 69)- (2. 79) to obtain 
[U~, A, r^] ♦ [\, D~, a, K] (5.43) 

Let 

X = K-G (5.44) 

Compute the updated factors [U~, D~] by applying the 
Agee-Turner triangular factorization (Appendix C) to 
the identity 


. 0jDi5<^T . 


aW 


(5-45) 


-I 

4 4 

I i 

i » 

< 


I I 


Covariance 

Time 

Update 


^xx* 


Let 


• «Stx - «xJ = ♦a -♦ 



""-Ba 

•» 

G 

[X, 

^XX 


_®a 

0 

1 

t U 

X 


D = diag (Qj^, Dj-, Djj) 


(5.46) 


(5.47) 

(5.48) 


Apply MWGS algorithm,^ Eqs. (3.28)-(3.30), to the 
W and D arrays to obtain 
[W, D] ♦ [U^f, UjQj, U,,, Djf, D^3 (5.49) 

At any stage, the true estimate error covariance is given by 


1 . I 


Py . . Uj^D,4x 


(5.50) 


^The modified Givens update, Eqs. (3. 80)- (3* 89 ) t oould also be used 
here. 
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This general error analysis technique inherits the computational 
efficiency of the 0-D filtering algorithms upon which it is based. 

The evaluation scheme is easily Implemented, requiring only minor modifi- 
cations in the U-D filtering equations. The following sections illustrate 
how this method is easily adapted to handle some important reductions 
of the general error analysis problem. 


5.3.2 Evaluation of Mlsmodeled Biases and Process Noise Parameters 


Consider the following application of the U»D error analysis 
technique. Suppose the parameter vector, X, is partitioned such that 


X . 


*2 ^"2 


(5.51) 


and let AA , A* , and * have the form 


"l "2 


AA s [0 I AA,] )1 


(5.52) 


n, n, 


A# a [0 i A^2] }n^+n2 


(5.53) 


"2 


1 *vA 


*a * — 


(5.5M) 
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For this problem the error equations corresponding to (5.23) and (5.25) 
reduce to the following. 


AX s (I-GA)aX + GAA2X2 + Gvg 


AX = *AX - A #2X2 - 


Hence, estimate errors contain only the X 2 portion of the actual 
state. Notice also that X 2 is uncoupled from by Eq. (5.54). Thus 
the error analysis corresponding to Eqs. (5.29)-(5.40) simplifies and 
involves only the vectors AX and X 2 . Let be partitioned compatibly 
with In this case the error equations (5.30) and (5.38) take the 


1 n^-«-n2 


n.|-i-n 2 { aX 


noC X 


— va }1 

G (I-GA)U5f Ujjj — 

Vv )n.+no 

0 0 u„ 

* J 


‘’x '"2 


k n.|+n2 H2 


n^+n2{ AX -Bg *U2 Ujjj -- 

n^{ X2 B 2 0 


}n^+n2 


where 


Ujfjj = (I-GA)U5f3j + GAAgU, 
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Thus, in Eqs. (5. 42)- (5. 47), the matrices A# and AA reduce to 
A*2 and AA2» respectively. The array in the lower portion of W, 
Eq. (5.47), reduces to B2 (cf 5.58). 


t ■ 

i' 


Suppose that X2 represents biases or colored noise parameters, 
so that *2 diagonal. In this case, each mismodeled component of 
X2 introduces only one additional "true" state into the error analysis; 
i.e., one additional row and column onto the augmented U and D arrays 
in Eqs. (5.42)-(5.49). Thus, filter sensitivity to mismodeled biases 
and colored noise parameters may be easily and efficiently evaluated. 


i 


u 

ia 

IS 

i 


5.3.3 Sensitivity Analysis of Onestimated Process Noise Parameters 
Another important reduction of the general error analysis 
problem Involves the evaluation of unmodeled colored noise parameters. 

In this case the estimate error recursions (5.29)-(5.40) are written 
in terms of Ax and p, where x represents the parameters estimated by 
the filter and p denotes the unmodeled parameters to be evaluated. 

The ♦gtAA and A» arrays in Eqs. (5.42)-(5.47) for this case reduce 
to 4A=Ap and A*s*jjp. Further modifications may be required 

in the W and D arrays (Eqs. 5.47 and 5.48), depending ^ the white 
process noise model used in the evaluation. 

The effects of unmodeled bias parameters could also be evaluated 
by this reduction of the general error analysis algorithm. However, 
bias effects are more efficiently evaluated by the sensitivity analysis 
technique described in the next section. 
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5.4 Variable Dimension Filtering and Sensitivity Analysis 


5.4.1 


f Unestlmated Bias Parameters 


When dealing with higher dimensioned systems involving 
large numbers of hl«>3-type parameters, it Is often necessary or desirable 
to neglect certain parameters in the filter model. Reasons for studying 
a reduced order problem include practical considerations of computer 
time and storage and concern that high-order filters can be overly 
sensitive to numeric effects. The following error analysis technique 
is of prime importance in the design and evaluation of reduced-order 
filters. 


Suppose the parcuneter estimation problem involves the vector 


X }n 


y }b 


The parameters x are to be estimated and suppose y are either a) esti- 
mated, b) neglected in the filter model, or c) neglected in the filter 
model but considered when the true covariance of x is computed. Let 

A 

the optimum estimate of the complete problem be X with error covariance 
where 0 and d are partitioned consistently with 


^ I L ^ ” 


Uy lib 


n b 

and s diag (D^, Dy) 


(5.61) 
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In terms of second-order statistics the 0-D factorization of P implies 
that estimate errors can be written as 


n{ r Ax 


b{ Ay 


«x «xy ''x 


(5.62) 


Vy )b 


The uncorrelated random vectors and have covariances anA Dy 
respectively. Equation (5«62) implies that 


a5 . 0,., . 


(5.63) 


Let ^ represent the estimate of x that corresponds to a n-dlmenslonal 
filter, i.e., A^ s 0 with probability one. If denotes the error 
covariance of then Eq. (5.63) Implies that 




(5.64) 


Note further from Eq. (5.63) that the **sensitlvlty" matrix, z, which 
relates the x estimates to errors in the y parameters^ is given by 


I = U^yUy* 


(5.65) 


Let ^ and ^ represent the optimal parameter estimates based on the 
complete model and let and ^y denote their respective covariances. 
Denote as the *'con.«'lder" covariance of i.e., the error 

covariance based on the complete model. Then Eqs. (5. 63)- (5. 65) yield 
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where 


^ + r ^ iT 

X y 

Moon) ^ ^ IP (o)jT 

* X y 


( 5 . 66 ) 

(5.67) 

( 5 . 68 ) 


and Py(o) Is the a priori covariance of the y parameters. The coveu*iance 
relations in Eqs. (5. 64)- (5. 68) were first recognized by Bierman [1976a]. 
He showed that ^ and are related as follows. 




(5.69) 


The reader may wish to verify from Eqs. (2. 69)- (2. 79) that 
and are computed independently and are not affected by the last 
m columns in the U and D arrays. 


This "consider" error analysis is easily extended to the 
discrete linear filtering problem studied in Chapter III if the transi 
tion equation has the form 



(5.70) 


where E{wi} s 0 and ECw^w*^} * diag (q*, ***iqi.). 
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The results of Chapter III imply that the U-D time update for this 
system may be accomplished by a square-root-free Givens triangularization 
of the array 


W s 


xy 


'y J 


}n 

}b 


(5.71) 


k n b 

involving the diagonal D s (Q, Dy). The elements of W in Eq. (5.71) 
are given by Wy = *yOy, and W^^y = + *xyOy 

Bqs. (3.75)- (3. 77) it is clear that after the last b columns of W and D 
are packed by the modified Givens algorithm we are left with 


,(n) 


xy 


'y J 


}n 

}b 


p(n) 


k n 


S (Q, D^, Dy) 


(5.72) 


(5.73) 


The and 5^^ factors are then obtained by applying the modified Givens 
algorithm further to the arrays W' s [B and D' = (Q, Dj^). Thus, 

and are computed in a way whioh is independent of the y parameters. 
For this reason the "consider** error analysis, Bqs. (5-6iO-(5.69) , 
applies directly to systems of the form (5.70). A separate analysis 
is required to evaluate the effects of neglected stochastic parameters 
such as colored process noise (cf section 5.2.3). 
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Bierman [1976a] has noted the following features of the 
"consider " analysis technique. 


1) From the results of a single large filter which includes 
all of the relevant parameters it is possible to 
evaluate a variety of filter models. Thus, beginning 
with the complete solution one can recursively compute 
estimates for models with n^^ parameters from estimates 
corresponding to the nj^^^<-dimenslonal model where 

”i ^ "i+1* 


2) Parameters may be rearranged for additional flexibility. 
That is, one may decide to estimate the last b param- 
eters and consider the first n. In that case the 
error equation (5.62) must be rearranged and has 
the form 



(5.74) 


This expression may be trlangularized by applying 
one of the U-D factorization algorithms in Chapter 
III. For example, a MWGS factorization may be applied 
to the arrays 


W ■ 



0 


and U « dlag (Dy, D^) 
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to yield 



(5.75) 


where v' and v' are uncorrelated and have covariances 
Dy and respectively. 

3) When the a priori covariance Py(o) in Eq. (5.67) 

is diagonal, the columns of £ P,,(o)^^^ represent 1c 
pe. .urbations in the estimate errors due to the corre> 
spending y parameter. This assoolatioi is useful 
in determining filter sensitivity to various parameters 
which may be omitted from the final filter model. 


4) The matrix inversion in Eq. (5.65) Involves a triangular 
matrix and, thus, can be easily accomplished by back 
substltut. on methods. 


5) The U-D variable-dimension filtering: method is more 

flexible, more reliable and considerably more efficient 
than corresponding techniques employing partitioning 
of the covariance matrices (of Frledland [1969]). 


5.4.2 gaag^dac, £UUc 

The consider filter to be discussed here is a generalisation 
of the method described by Schmidt [1968]. The filter model includes 
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parameters x and y idiere only the x variables are estimated. The y 
parameters and their associated error covariance are retained at the 
a priori values. This method differs from the consider error analysis 
studied In the previous section since the gain, G^, used to update 
the X parameters Is computed from the full "true** covariance on x and y. 

This type of suboptlmal filter might be used, for example, 
in situations id;ere the y parametera cannot be accurately estimated 
or when estimating them might lead to an overly optimistic y covariance. 

In the latter situation the filter is, in effect, reduced to an x filter. 
Mhen It Is believed that a reduced-order filter cannot accurately estimate 
X, the consider filter is sometimes proposed (cf Bierman [1976b]). 

Error propagation for the consider filter takes the form 



where 

z * AjjX ♦ kyY * 
y > B{y) « 0 


(5.76) 

(5.77) 

(5.78) 

(5.79) 
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If the U-D factorization is assumed for the true covariances P and 
P associated vilth these errors, then Eq. (5.76) may be Interpreted 
as follows 


n{ 

b{ 


n 


b n b 


0 





( 5 . 80 ) 


where v^, and are unoorrelated with diagonal covariances Dy 
and r^, respectively. Equation (5.80) may be rearranged to yield 
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( 5 . 81 ) 


( 5 . 82 ) 


The covariance update corresponding to (5.76) may be Interpreted as 
a triangularizatlon of the error expression (5.81) such that 
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where 0^^ and are uncorrelated with diagonal covariances 6^^ and 6y 
respectively. The constraint that 


^ s B' 






I 


be unchanged by the transformation from (5.81) to (5.83) yields 


[Oy, fty] = [Oy, Dy] (5.84) 

®xy * °xy ^?*®5) 

= <I-6xAx>6x5xOx^(I-GxAx)^ ♦ GxVx" <5.86) 


If were the Kalman gain associated with the r'^duccd* 
order problem having a priori covariance factors Uj^ and 5^^, then the 
factors Ojj and could be obtained by applying the Bierman 0-0 measure- 
ment update algorithm to Eq. (5.56). Note, however, that for the general- 
ized Sohmldt consider filter, the gain 0^ Is computed from the full 
consider covariance. Thus, la obtained from the identity 



a 



(5.87) 


( 5 . 88 ) 
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fX> 


(5.89) 


(5.90) 


For this reason, we compute 0,^ and by applying the U-D arbitrary 
gain update algorithm (section 5.2) to Eq. (5.36). The first phase 
of this update applies the Bierman algorithm, Eqs. (2.69)-(2.79) , to 
obtain the gain K = corresponding to the reduced-order problem 

with a priori factors Ujj and 5,^. Note that the full-dimension gain 
in Eq. (5.87) is given by G = Kn+b^“n+b ^n+b “n+b **® 

obtained by cycling further through Eqs. (2.7^0 and (2.78) for 
j = n+1,...,n+b. In this recursion only the first n elements of the 
vector Kj need be computed since the gain Gy is not required. 


Time propagation via (5.77) of the full consider covariance 
U-D array is easily accomplished by a MWGS or modified Givens factoriza- 
tion applied to 
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Chapter VI. Mumerlcal Comparison of Discrete 
Kalman. Filtering Algorithms 


6. 1 Introduction 

The numerical stability and accuracy of various discrete Kalman 
filtering algorithms have been carefully studied by applying these 
methods to a realistic interplanetary navigation problem. Included 
in the study were the conventional and stabilized Kalman algorithms, 
the Potter-Schmidt square root filter, and the U-D factorization filter. 
The Carlson-Schmldt algorithm was omitted from the comparison since 
this method clearly shares the numerical characteristics of the U-D 
filter (ef Chapters II and III), 

Each of the algorithms was implemented on the UNIVAC 1108 computer 
vrtiich has a 60 bit characteristic (18 decimal digits) in double precision 
and 27 bit characteristic (8-9 decimal digits) in single precision. 

The complexity of the study problem prohibits closed form solutions. 
However, numerical solutions from all the algorithms, using double 
precision arithmetic, agreed to 8 digits or more. These double precision 
results were used as a reference for evaluating the accuracies of the 
same algorithms computing in single precision. 

6.2 Pr.2fcl.gffl.Dg.9firiBtl9r> 

A portion of the 1977 Mariner Jupiter-Saturn (MJS) deep space 
mission was chosen for the filter comparison study. The problem involves 
spacecraft navigation for the 30 day period immediately preceding 
Saturn encounter (point of closest approach). For the initial 20 days 
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the spacecraft trajectory is nearly rectilinear, a characteristic which 

is typical of a major portion of most deep space missions. The last 

segment of the trajectory has a hyperbolic bend due to the gravitational 

effects of Saturn. Hence, tracking data from this phase of the mission 

is particularly useful for an accurate determination of the planetary 

mass and the spacecraft position and velocity. This trajectory is thus 

characteristic of a significant class of, orbit determination problems. 

* 

The nominal trajectory and state transition matrices were obtained 
by integrating the equations of motion and variational equations (cf 
Moyer [1971]). Earth-based measurements of the spacecraft Include 
both doppler and ranging data. Partial derivatives of the data with 
respect to the system parameters were evaluated analytically about 
the nominal trajectory using the orbit determination software described 
by Moyer. All of these calculations were performed in an earth-centered, 
cartesian coordinate system. 

Perturbations from the nominal trajectory and simulated data 
were constructed from a discrete linear model. The complete model 
Includes differential corrections in 19 parameters; the spacecraft 
position, velocity and acceleration components (3 each); the gravitational 
constant of Saturn (GM); and tracking station locations (3 cartesian 
components for each of 3 stations). The state transition for these 
perturbations is described by the following equation. 
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( 6 . 1 ) 


The vector x denotes the 6 components of position and velocity, p repre- 
sents the 3 acceleration components, and y represents the 10 bias parame- 
ters. Accelerations are modeled as piece-wise constant colored noise 

—11 2 

with time constants of 12 hours and standard deviations of 10 km/sec ; 
and these define the variances of the white noise, w^, in Eq. (6.1). 

The discrete times (tj^} are assumed to be equally spaced with 

= At taken as 2 hours. Thus has covariance Q = ql^ where 

q = (1 - m^) Op^ = (1 - e"^'^^) 10"^^ (km/sec^) (6,2) 

Let X represent the complete state vector in (6.1). T'.te initial 
vector Xp has zero mean and diagonal covariance Pp. The nontrivial 
elements of Pp are defined by the standard deviations in Table 6.1. 

These a priori statistics are typical assumptions for the kind of orbit 
determination problem considered here (of 0"Neil, et al. [19733) • 

The simulated state Xp is obtained from a Gaussian random variable 
generator with zero mean and covariance Pp.*^ The actual state is then 
propagated according to Eq. (6.1) where the components of w^^ are obtained 
from the Gaussian random variable generator with covariance Q. 

^Eaeh component of Xp is obtained by using a N(0,1) ranaom number gen- 
erator followed by an appropriate scaling of the sampled value. 
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Table 6.1. A Priori Statistics Used to Generate Sample Path 



Data for the linear simulation are generated from the equation 


H = *i*i + ''l 


where the elements of are partial derivative coefficients evaluated 
along the nominal trajectory, and is white data noise obtained from 
a Gaussian random number generator using the appropriate measurement 
error covariance. The simulation includes two or three doppler measure 
ments every two hours with 1 nim/sec accuracy (for 1 minute averaging 
time) and occasional range points with an accuracy of 3 meters. There 
were a total of 535 doppler and 72 range measurements in the 30 day 
tracking period. 
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Thlj simulation is representative of a large number of interplan- 
etary navigation problems. The a priori statistics are in no way 
contrived to produce poor results in the conventional algorithms. On 
the contrary, a priori state covariances were chosen on the small side 
for problems of this type in order to avoid initialization errors asso- 
ciated with the Kalman algorithms. Similarly, process noise uncer- 
tainties were sissumed to be an order of magnitude higher than usual 
for this kind of mission^ because the Kalm 2 ui algorithms generally expe- 
rience less numerical difficulties in high process noise environments. 

This estimation problem is well posed in an engineering sense. 

The problem is observable; the transition matrix is not ill-conditioned; 
and the measurement coefficients and a priori error variances are not 
unusually large. Thus, the results of this study should be of interest 
to the entire estimation and control community. 

6.3 Filter Implementations 

The five covariance-type filter algorithms compared in this study 
are the conventional Kalman filter, Joseph's stabilized Kalman filter, 
a conventional Kalman filter with lower bounding, the Potter-Sohmidt 
square root filter, and the U-D factorization filter with MWGS time 
propagation. A limited evaluation of the U-D filter using the modified 
Givens update is also Included. Details of the various filter algorithms 
are given in Chapters II and III. 


^See Christensen [1976] or O'Neil, et al. [19731. 
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Each algorithm Involved in the study was mechanized for maximum 
accuracy and efficiency. Thus, the U-D filter employs the efficient 
one-at-a-time colored noise update, Eqs. (3.130)-(3.133); while the 
Potter-Schmidt filter performs time updating in one stage, taking ad- 
vantage of special system structure and blocks of zeroes. 

Computations are also reduced by using vector outer products 
whenever possible. In the Kalman algorithms, symmetry of the covariance 
matrix is preserved by computing only the upper triangular elements. 

The stabilized Kalman filter contains a single exception to this rule. 

Significantly improved results were obtained with the stabilized 
Kalman filter when off-diagonal covariance elements were computed and 
averaged during a portion of the measurement updating. This averaging 
takes place during the computation of the first term in Eq. (2.24). 

Thus all n^ elements of the array P = - V 2 K*^ are computed, and 

the off-diagonal elements of f are obtained as follows. 

^ij = + P ,' + rKj^Kj i ^ j (6.4) 

The fact that numerical re^ are sensitive to such mechanization 
details is indicative of the algoi .uim's instability. Even with the 
averaging of off-diagonal elements, the stabilized Kalman algorithm 
performed poorly. 

The conventional Kalman filter with lower bounding is designed 
to perform a covariance measurement update in the following way. The 
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conventional update formula is first applied to obtain an intermediate 
array P. Thus, 


P = P - KAP 


(6.5) 


The updated covariance P is then computed from Eqs. (6.6)-(6.8). 


*JJ “ ^ " I*****'' 


( 6 . 6 ) 


If fu ^ "ij 


sgn(Pj^j) otherwise 


(6.7) 


where 


2 A A 
^ij = Pmin **ii 


1 = 1 , 2 ,. 


( 6 . 8 ) 


The n components of and the correlation are chosen 
a priori. This mechanization is typical of the techniques that are 
frequently used to prevent the computed covariance from becoming lndefl« 
nlte (of Schmidt et al. [1968]). The mechanization is no optimal 
and the computed ^ is generally not the actual error covariance. Choosing 
the bounds and la something of an art, and appropriate values 
are usually determined from lengthy simulation studies. This lower 
bound filter algorithm is included in the comparison study in order 
to illustrate that ad hoc "patching" techniques can compensate to some 
extent for the numerical inadequacies of the conventional algorithm. 
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However, the results of the remaining sections demonstrate that such 
measures are not necessary when factorization algorithms are used. 

All of the filter algorithms were coded and carefully checked 
using double precision arithmetic. These programs were then converted 
to single precision by removing the FORTPAN IV "implicit double preoisioa" 
statement. However, in the single precision programs, inner products 
are accumulated in double precision and then rounded to single precision. 
In addition, estimates are retained in double precision because a signifi- 
cant accuracy loss is Incurred when single precision arithmetic is 
used in the estimate propagation, Eqs. (3.11) and (3* 13). This accuracy 
loss would not be attributable to a particular filtering method but 
would mask the desired comparison. 

Since all of the filter programs compute estimates in double 
precision and use the same (single precision) inputs,^ only the gain 
computations differ. Thus, filtering accuracies may be analyzed by 
applying the gain evaluation algorithm described in Section 5.2. This 
evaluation program, implemented in double precision and operating with 
the same inputs used by the filter, computes the actual, or true, error 
covariance associated with each gain profile; i.e., each filter algorithm. 
Hence, two complimentary methods are used to compare filter accuracies. 

The error analysis program provides a statistical means of comparing 


’^The state transition and data matrices used in the simulation are 
computed in double precision and rounded to single precision. 
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filter performance, while the single simulation results yield an illustra- 
tion of these statistical predictions. 

6.4 Numerical Results 

A comparison of filtering accuracies was first obtained for the 
full 19-state model described in section 6.2. The filter model was 
then varied In order to evaluate sensitivities to a priori statistics. 
Reduced-dimension problems were also studied as a means of assessing 
how numerical stability is affected by model dimensionality. In every 
case studied, the double precision programs computed estimates and 
standard deviations that agreed to 8 digits or more. In fact, the 
factorization algorithms agreed to at least 10 digits. The single 
precision programs, however, produced widely varying estimates and 
covariances. The results to be described indicate that, while the 
single precision Kalman algorithms experience serious numerical deteriora- 
tion as computer word length decreases, the factorization methods have 
excellent precision and stability. 

6.4.1 The Complete State Model 

For this portion of the study, all of the filters assumed 
the 19-state model described In section 6.2. Thus, the actual and 
assumed models coincide. For this case and the others to follow, the 
single precision factorization algorithms computed gains and variances 
which agreed with the double precision reference results to about 5-6 
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digits . ^ The single precision Kalman algorithms, however, produced 
markedly different gain profiles and variances. Numerical deterioration 
of these algorithms was evident In the form of negative computed variances 
which appeared at Inexplicable times. Although this numerical Instabllty 
Is related to the choice of a priori statistics (cf section 6.4.2) 
negative variances do not appear until after four days cf filtering 
with the conventional algorithm, and after ten days when the stabilized 
formula Is used. This behavior Is particularly Interesting when we 
note that In the first four days there are 48 time and 80 measurement 
updates. 

Both the conventional and the stabilized algorithms compute Intermit- 
tent negative variances. Even the bias parameter variances are Intermit- 
tently negative. Furthermore, these negative elem'ints are not related 
to filter variances which are approaching zero. The erratic behavior 
of the stabilized algorithm Illustrates this point. For example, at 
9.75 days the stabilized formula computes « -1.8 x 10^ (km^/sec^)^ 
and at ten days adjusts this to 1.7 x 1o\ The correct (double precision) 
value, however, is ~ 5 x 10^ (kmVsec^)^. 

The extent of the numerical deterioration In the Kalman algorithms 
Is apparent In Figures 6.1 and 6.2. These figures compare the actual 

^The U-D filtering results reported In this section were obtained with 
both time update methods; the MWGS and the modified Givens algorithms 
(of sections 3.3 and 3.4). The latter propagation scheme employed 
Method B exclusively (of Equ. (3.80)-(3.89)) and experienced no numerical 
degradation; l.e., both filter Implementations produced results which 
agreed to about 5-6 digits with the double precision standard. 
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position and velocity uncertainties^ associated with each filter and 
were obtained by applying the gain evaluation program. 

In Figures 6.1 and 6.2 the U-D and Potter single precision filtering 
accuracies are shown to agree with the dc’)hi« precision references. 

Notice » however, the excursions of the conventional and stabilized 
Kalman single precision accuracy curves. The conventional Kalman results 
are seen to be particularly inaccurate at 4 days, while the stabilized 
algorithm deteriorates rapidly at 10 days. 

Simulated data, generated using the same model assumed in the 
error analysis, was filtered by each of the algorithms. Actual estlmr:<&.'‘ 
errors for this single simulation are compared in Figures 6.3 and 6.4. 
Notice how closely the gain evaluation results in Figures 6.1 and 6.2 
predict the error curves of the sampls path. As anticipated, the conven- 
tional algorithm in single precision produces large errors at 4 days. 

The single precision stabilized algorithm also performs poorly, especially 
after 10 days. However, the estimates and computed standard deviations 
obtained from the stabilized filter, when monitored at one day Intervals, 
show few signs of numeric deterioration. Except for the times (3) 
that negative variances are printed, these estimates and statistics 
appear reasonable and oonsistent. Only when the results are oompared 


^The root-sum-square of position and veloolty errors are chosen as 
a measure of estimation accuracy beoaune these parameters are of pri- 
mary Interest in navigation and are representative of the general 
filtering results obtained in this study. 
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with the double precision references is it apparent that the stabilized 
Kalman estimates are far from optimum. 

Recall that the stabilized formula was introduced as a computational 
improvement to the conventional method. It is, therefore, surprising 
to note that after 10 days the conventional algorithm yields significantly 
more accurate results than the stabilized formula. This phenomenon 
was not observed in all of the cases studied, however. The point to 
be made here is that, although the stabilized algorithm does generally 
give improved performance over the conventional formula, neither method 
can bj considered reliable. 

By comparison, estimates computed using the factorization algorithms 
agree to about 4 or 5 digits with the double precision values. This 
agreement corresponds to better than 1 km in position and 50 mm/sec in 
velocity. These single precision accuracies are particularly impressive 
vrtien it is noted that estimation uncertainties are two orders of magnitude 
greater than these differences (cf Figures 6.1 and 6.2). In other 
words, the differences in the single and double precision results are 
in the "noise" level. 

Jn every case studied, the relative position and velocity accuracies 
displayed the same general agreement illustrated in Figures 6. 1-6.4. 

In view of this, subsequent discussions are restricted to the comparison 
of position errors and uncertainties. Thus, unnecessary discussions 
of velocity results are omitted. For similar reasons, the conventional 
Kalman algorithm is also omitted from subsequent discussions. The 
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numerical instability of the conventional algorithm is already well 
documented (cf Bellantoni and Dodge [1967] or Schmidt [1968]). This 
method is dismissed from further comparisons with the observation that 
in nearly every case studied the single precision conventional formula 
suffered numerical degradations similar to those displayed in Figures 

f 

6. 1-6.4. 


Numerical divergence of the Kalman filter is generally associated 
with Indefinite computed covariance matrices. Hence, it is common 
practice to attempt to preserve positivity by bounding the diagonal 
elements from below and to limit the correlations between pairs of 
variables (cf Bqs. (6.6)-(6.8)) . Any such attempt to stabilize the 
conventional Kalman algorithm Introduces a myriad of filtering alterna- 
tives. For example, should the lower bound on the velocity uncertainties 
be 1.0 m/sec or 0.1 m/sec, and should the maximum correlation be .99 
or .98? The choice of such patch factors is often problem dependent 
and may require lengthy simulations. 

In this study filtering results are Indeed sensitive to the choice 
of bounds as Figure 6.5 illustrates. This figure displays the RSS 
position error profiles produced by the single precision patched algorithm 
for various bounding schemes. By comparing Figures 6.3 and 6.5 one 
can conclude that patching yields a marked Improvement over the stabilized 
Kalman results. However, all of the error curves associated with the 
patched algorithm are far above the optimal double precisinn results 
in Figure 6.3> Actual filtering uncertainties a>e compared in 
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Figure 6.6, and the patched algorithm is shovm to be significantly 
less accurate than the factorization filters. 

In all of the cases studied, the patched Kalman filter displayed 
the same poor performance observed in Figures 6.5 and 6.6. This study 
demonstrates that the practice of introducing ad hoc patch factors 
to combat Kalman filter numerical divergence results in algorithms 
that are significantly less efficient and accurate than the factorization 
methods. Patching techniques are thus omitted from further consideration. 

6.4.2 Effects of Incorrect A Priori Statistics 

Numerical difficulties with the Kalman algorithm can often 
be attributed to initial ill-conditioning caused by large a priori 
state uncertainties and relatively small data covariances. These problems 
can be reduced by scaling the a priori statistics. However, any improve- 
ment in numerical conditioning is offset somewhat by the effects of 
suboptimal modeling. Consider, for example, the case where Initial 
velocity uncertainties are reduced by an order of magnitude to 10 m/sec, 
and range uncertainty is increased to 10 meters (from 3 meters). This 
combination of a priori statistics, selected by numerical experimentation, 
allows the stabilized Kalman algorithm to appear "stable." In fact, 
for this choice of filter statistics neither the conventional nor the 
stabilized algorithm computes negative variances. Moreover, for this 
example simulation estimate errors are consistent with the filter formal 
statistics. This consistency creates the false impression that the 
Kalman algorithms are performing well when, in fact, they are grossly 
Inaccurate. Actual Kalman accuracies, computed by the gain evaluation 
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Pig. 6.5. Comparison of Patched Kalman 
Filter Errors for Different Reset Values, 
19~State Single Simulation Results 



Fig. 6.6. Performance Comparison of 
Patched and Stabilized Kalman Filters, 
Complete 19~State Model 
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program using the correct statistics, are displryed in Figure 6.7. 

This figure shows that for most of the filtering period the stabilized 
Kalman uncertainties (using scaled a priori statistics) are an order 
of magnitude greater than achievable filter accuracies. 


The results in Figure 6.7 also show that the Kalman filter is ' 
more accurate when suboptimal (a^ = 10 m/sec, og = 18 meters) rather ■ 
than optimal (oy = 100 m/sec, oj) > 3 meters) statistics are assumed. > 
This theoretical impossibility is, of course, due to numerical errors i 
which can cause unpredictable results. i 


1? 

4 
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k 
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When only one of the a priori uncertainties (oy or op) is scaled, 
the single precision Kalman filters continue to produce negative computed 
variances and unreliable gain profiles. Actual filtering accuracies 
corresponding to scaled a priori velocity variances are illustrated in 
Figure 6.8. In this case, Oy is scaled down by an order of magnitude 
in the filter model, and gain profiles are evaluated using the correct 
statistical model. Note that, for this example. Initial velocity vari- 
ances are scaled down by two orders of magnitude. However, instead 
of Improving Kalman filter accuracies, this scaling produced greater 
errors. Figure 6.8 demonstrates that the stabilized filter (with 
Oy s 10 m/sec) yields large uncertainties at 5 days which peak again 
at 22 days. By comparing Figures 6.1 and 6.8 we can conclude that 
numerical errors in the Kalman results occur at completely arbitrary 
times and are unrelated to any physical phenomena peculiar to the 
problem. 
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Fig. 6 . 7 . Performance Comparison 
of Stabilized Kalman Filters Using 
Scaled A Priori Statistics, 19-State 
Filter Model 



Fig. 6.8. Comparison of Actual 
Position Uncertainties for Scaled 
Velocity A Priori, 19-State 
Filter Model 



Fig. 6 . 9 . Comparison of Position Errors for Scaled 
Velocity A Priori, 19-State Single Simulation Results 
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In a filtering problem with observability and significant amounts 
of process noise, one would expect estimates to depend primarily on 
data in the recent past. Thus, estimate error profiles corresponding 
to different a priori covariance assumptions should, except for initial 
transient effects, look quite similar. The factorization filters illus- 
trate this effect in the bottom two curves of Figure 6-8. In every 
case involving scaled a priori statistics, the factorization filters 
demonstrate the same stability evident in Figure 6.8. The Kalman algo- 
rithm, or: the other hand, is quite sensitive to the choice of a priori 
statistics as the topmost curves in Figures 6.7 and 6.8 illustrate. 

The single simulation results for this case yield estimate 'ror 
profiles close to those predicted by the error analysis. Position errors 
are compared in Figure 6.9. Notice the striking resemblance between 
the factorization filter error curves in Figures 6.3 and 6.9, particularly 
after 10 days. As expected from the error analysis, the Kalman results 
in these two figures appear totally unrelated. 

This analysis illustrates how numerical instability can cause 
unpredictable results which violate established estimation principles. 

6.H.3 Reduced-Dimension Problems 

The discussions in the previous sections were based upon 
the complete 19-state model. In this section models of smaller dimension 
are examined. This study shows, among other things, that the numerical 
instability of the Kalman algorithms is not caused by the dimensionality 
of the model; and that the inclusion of process noise improves the 
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appearance of the computed covariances but not the accuracy of the 
estimates. 

The smallest, physically meaningful model corresponding to the 
planetary approach problem Involves the six position and velocity vari- 
ables. This six-state system represents a parameter estimation problem 
since there is no process noise. The Kalman algorithms are known to 
be numerically unstable for parameter estimation problems. Hence, 
it is not surprising that the single precision Kalman filters compute 
intermittent negative variances when the six dimensional model is assumed. 
However, these 6-state filters, when applied to the simulated data, 
manage to at least partially track the spacecraft. Actual filtering 
accuracies for this case are displayed in Figure 6.10. The top two 
curves in this figure represent the accuracies obtained with the six 
dimensional filters operating in the complete 19-state environment. 

Notice that the stabilized Kalman algorithm suffers a severe accuracy 
degradation in this case. Position uncertainties are two orders of 
magnitude greater than those obtained with the correct model (cf Figure 
6.1). By comparing the factorization results in Figure 6.10 (bottom 
two curves) one can see that the accura'*y loss due to m.ismodeling is 
considerable . The stabilized Kalman results in this figure suggest 
that either this algorithm compounds the effects of mlsmodellng or 
numerics represent the dominant error source. 

Numerical effects may be further separated from mode, ag errors 
by evaluating filtering accuracies for an actual slx-dim«»'' jnal model 
Hence . no effects due to mlsmodellng are present in the .uracy 
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comparisons of Figure 6.11. Notice that position uncertainties corre- 


sponding to the stabilized algorithm are nearly Identical in Figures | 

6.10 and 6.11. These figures demonstrate that numerical errors asso- | 

j 

dated with the stabilized algorithm are so large that they completely ’ 

obscure the effects of mismodeling. By contrast, the factorization j 

curves in Figures 6.10 and 6.11 clearly indicate how 6-state filtering 
accuracies are degraded by the presence of unmodeled parameters. 

i 

\ 

The second model selected for careful study included the three 
colored noise accelerations in addition to the six position and velocity 

! 

parameters. This 9>state system was chosen because it Includes a slgnifl- ; 

cant amount of process noise, and it is generally assumed that high 

I 

levels of process noise will stabilize the Ka^^an filter computations. , 

The stabilized Kalman results appeared to support this theory. 

That is, the stabilized algorithm computed covariances, gains, and 
estimates (based on the simulation sample) which looked reasonable. 

No hint of numerical difficulty was evident. The results, however, 
differed from those obtained with the U-D and Potter-Schmidt filters. 

As usual the factorisation algorithms agreed closely with the double 
precision reference results. 


Error analysis for this case produced results similar to those 

i 

observed in the six-dimensional problem. That is, numerical errors 

in the Kalman calculations were again substantial and completely obscured I 

any effects due to unmodeled parameters. Figure 6.12 shows how severely : 

the Kalman results are degraded by numerical errors. This figure displays \ 
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actual position accuracies associated with the 9-state filters operating 
In an actual 9-state environment. These results demonstrate that numer- 
ical deterioration In the Kalman filter translates Into position errors 
which are an order of magnitude higher than they need be. 

This example Illustrates that, while the Inclusion of process 
noise Improves the performance of the Kalman algorithms (cf Figures 

i 

f 

I 6.11 and 6.12), they still lack the accuracy achievable with the factorl- 

f 

I zatlon methods. 

f 

i 

6.5 Conclusions 

This study demonstrates the excellent numerical stability and 
precision of the U-D and Potter-Schmldt factorization algorithms. 

Both methods. Implemented In single precision, produced results which 
were close to the double precision references. The numerical stability 
of these algorithms was further demonstrated by their lack of sensltlv.lty 
to the choice of a priori variances and process noise levels. 

i 

I The Kalman filters, on the other hand, were acutely affected 

by the use of single precision arithmetic and scaled a priori statistics. 

; Wampler [1970] observes that these are sufficient reasons to declare 

an algorithm numerically unstable and to abandon It. The results of 
this study support his conclusion. Both the conventional and the stabi- 
lized algorithms experienced severe numerical deterioration In nearly 
every case examined. Covariance matrices with negative diagonal elements 
were a common occurrence, and gain profiles were often erratic and 

! 

Inaccurate. 
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Even when the assumed statistics v.ere modified to stabilize the 


calculations, the single precision Kalman algorithms performed poorly. 
Accuracy degradation was often not apparent but had to be detected 
with a double precision gain evaluation program. 


In every case studied the U-D and Potter-Schmidt algorithms out- 
performed the Kalman methods. Accuracy improvements were generally 
substantial, often reaching orders of magnitude. Error analysis showed 
numerical degradation to be the dominant source of error in the Kalman 
algorithms. In fact, numerical errors completely obscured the effects 
of mismodeling. This result is of special interest because numeric 
effects are rarely considered when mission design requirements are 
constructed. These effects appear to be inconsequential, however, 
when the factorization algorithms are employed. 


The cost comparisons of Chapter IV show that for most filtering 
problems the factorization algorithms are not unreasonably costly. 

The Potter-Schmidt costs are generally within 80J of the conventional 
Kalman costs. The U-D algorithm, however, is competitive with the 
conventional method, and for some applications is somewhat more efficient. 
For problems involving a few colored noise parameters and a large number 
of bias parameters, the U-D algorithm is particularly efficient relative 
to the other methods. The complete 19-state model studied in section 
6.4.1 is an illustrative example. The CPU times for this case are 
compared in Table 6.2. These times Include the costs associated with 
Indexing, logic and other overhead costs not included in the analytic 
comparisons of Chapter IV. Table 6.2 shows the Potter-Schmidt algorithm 
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to be the most expensive filtering scheme, requiring considerably more 
time than the other methods. In sharp contrast, the U-D algorithm 
is seen to be the most efficient method, faster than even the conventional 
Kalman algorithm.^ 

This analysis of a meaningful engineering problem has demonstrated 
how the U-D amd Potter-Schmidt factorization algorithms can dramatically 
reduce the effects of numerical errors. Moreover, the U-D method combines 
numerical stability with an efficiency comparable to that of the origin . 
Kalman algorithm. 


Table 6.2. Comparison of Filter Execution Times* 
for the Complete 19-State Model 


Filter 

Single Precision 

Double Precision 

Conventional Kalman 

39 

n9 

Stabilized Kalman 

ns 

59 

U-D (MWGS) 

38 

46 

Potter-Schmidt 

63 

80 

*CPU time in seconds. 


These times do not include the costs of variance calculations in the 
factorization filters. However, even if variances were computed at 
every stage, the total CPU time would be increased by no more than 
ten percent. 
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Chapter VII. Sumnary of Contributions ^nd Areas 
for Future Research 

7.1 Summary of Contributions 

A new computational form of the discrete Kalman filter has been 
developed. This method is based upon the U-D factori 'ation of the error 
covariance matrix and has the characteristics of improved accuracy and 
stability generally associated with square root filters. The U-D filter 
was obtained by extending the Bierman [1976a] data processing method to 
allow for time propagation. Two numerically reliable propagation algo- 
rithms were derived in Chapter III. These algorithms were obtained from 
modifications of the Givens and Graiu-Schmldt matrix trlangularizatlon 
methods which are known for their numerical stability. Th"- general U-D 
propagation results were then extended to derive an efficient U-D time 
update algorithm for systems with bias parameters and colored process 
noise . 

Cost comparisons in Chapter IV have demonstrated that both cf the 
basic U-D propagation methods are efficient and usually require equi- 
valent amounts of computation. These comparisons further show that the 
U-D filter algorithm, using either the MWGS or the audlfied Givens 
propagation scheme, has efficiency nearly equal to that of Kalman's 
original formula. Moreover, for problems involving modest r umbers of 
process noise parameters, the U-D method is somewhat faster than the 
Kalman algorithm. By comparison, the square root covariance filters 
were shown to be significantly less efficient. In tyoical real-time 
applications Involving small computers and low-order systems, the square 


I 
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root algorithms are at least 40j more time consuming than the conven- 
tional Kalman method and may require in excess of 100$ more computation 
time. 


An extension of the U-D filtering results has produced a new co- 
variance error analysis technique. This method was derived in Chapter 
V as a mecins of evaluating general modeling errors, sue i as incorrect 
a priori statistics, incomplete parameter sets, and mismodeled system 
dynamics and data equations. The general error analysis algorithm re- 
quires only a modest alteration to the U-D filter equations. Efficiency 
is competitive with that of conventional error analysis methods, and 
numerical accuracy of the results is believed to be improved. The U-D 
evaluation method is easily adapted to analyze a variety of important 
special problems, such as the mismodeling of colored process noise and 
bias parameters. In addition, this method has produced a new consider 
f Iter algorithm which shares the simplicity and efficiency of the opti- 
mal U-D filter equations. 

The orbit determination case study in Chapter VI provides a thorough 
and extensive examination of the numerical characteristics of the various 
Kalman filtering methods. Ti._o study demonstrates with a meaningful engi- 
neering problem the improved precision of the U-D and square root covari- 
ance methods. -e methods consistently outperformed the conventional 
and stabilized Ki .m algorithms. Accuracy improvements were usually 
substanti;.! and often involved orders of magnitude. While the Kalman al- 
gorithms consistently experienced severe accuracy degradations and were 
generally unreliable, the covariance factorization methods exhib..tad 
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excellent numerical stability and precision throughout the study. A 
comparison of actual CPU times for this analysis has demonstrated that 
the U-D euid conventional Kalman filters require equivalent amounts of 
computation. Thus, this study illustrates how the U-D filter ccmbines 
superior numerical precision with exceptional computational efficieucy. 

7.2 Areas for Future Research 

The modified Givens algorithm which applies the efficient formula, 
Method B, has been recommended for U-D propagation. This algorithm is 
known to be numerically stable if the ratio is tested at each step 

and if the more stable formula. Method A, is used whenever d^/d_ is 
sufficiently large (cf section 3>^)< However, for most time updating 
applications of this algorithm, such testing is probably not necessary. 
Further analysis is needed to determine which problems, if any, are 
ac* ’ally susceptible to error growtl when Method B is used exclusively. 

Another reliable method for propagating U-D covariance factors can 
be obtained by appropriately modifying the Householder triangularization 
technique. Preliminary analysis indicates tnat square root calculations 
are Inherent to Householder techniques, and so this method may have 
little advantage over the MWGS or modified Givens propagation algorithms. 
However, a square- root- free Householder triangularization method could 
yield a U-D time update algorithm which requires fewer multiplications 
and divisions than the modified Givens formula. Hence this Householder 
method might be noticeably more efficient than either of the currently 
available U-D propagation techniques, particularly on computers where 
divide operations are relatively time consuming (cf section 4.3.1). 
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The U-D covariance error analysis method is based upon a number of 
algorithms which have been established as numerically reliable formulas. 
While the error analysis technique appears to in'.erit the stability of 
these individual algorithms, the integrity of this method should be con- 
firmed by further research. 

There is another area of covariance error analysis where research 
might prove fruitful. Although much effort has been expended to develop 
various gain evaluation techniques, very little attention has been given 
to the analysis of innovation computations. The case study in Chapter VI 
of this thesis has shown that parameter estimates are particularly sensi- 
tive to innovation errors. Recall that in the single precision U-D 
filter it was necessary to compute estimates and innovations in double 
precision arithmetic, although gain calculations were reasonably accurate 
in single precision. Further evaluation of Innovation errors could pro- 
duce some Interesting and useful results. 

Another extension of the U-D filtering method which would be 
useful, and might prove to be challenging, is the development of algo- 
rithms to allow for cross-correlations between measurement errors and 
process noise parameters. At first glance this cross-coupling appears 
to make the U-D time update unreasonably complicated. However, the 
"weighted error" propagation technique used in Chapter V could be applied 
to this problem and should make it apparent which orthogonal transforma- 
tions are required to accomplish the update. 
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The covariance factorization techniques presented in this report 
might e beneficially applied to problems involving continuous parameter 
systems and continuous data. These methods could also prove valuable in 
distributed parameter problems and in the areas of identification and 
control . 
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Appendix A 

Cholesky Decomposition and Data Whitening 

A. 1 Souare-Root-Free Choleskv Decomposition 

It has been shown (cf Martin et al. [1965]) that any positive 
definite symmetric matrix R can be uniquely factored as 

R s ODU*^ (A*1) 


tdiere U is unit upper triangular'*' and D is a positive diagonal. The U-D 
factors of R may be computed as follows. 


For j = n,...,1 evaluate recursively Eqs. (A.2)-(A.5). 


= r 


Jj 


n 

- L 

k=j+1 


Vjk 


(A. 2) 




(A. 3) 


"ji 


E 

k=j+1 


Vik'^Jk^l 


i = 


= 0 


(A.il) 


(A. 5) 


"^The factorization given here is actually a modification of the algorithm 
by Martin et al. [1965] since their formula computes a lower triangular 
factor. 
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A. 2 Choleakv Souare Root Decompoaltlon 

Let an upper triangular matrix S be defined aa 

S = ^ g) 

where ia the poaitive aquare root of 0 in Eq. (A.1). Equationa 
(A.1) and (A. 6) then imply that any poaitive definite aymnetric matrix 
R may be uniquely factored aa 

R = SS^ (A.7) 

> 

where S la upper triangular. Riatorically, thia decompoaition preceded 
the U-D algorithm and waa derived by Choleaky. The following algorithm 
ia a modeat rearrangement of the Choleaky aquare root factorization 
given in Fox [1954], 


For j = n,...,1 compute recuraively Eqa. (A.8)-(A. 10) . 


1/2 




Hu ®lk^ 
k=J+1 '' 


(A. 8) 



Table A. 2 contains a summary of the arithmetic operations required 
by this algorithm. 
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Table A. 2. 

Arithmetic Operations Required 


for a 

Square Root Decomposition 




Square i 

Computation 

Adds Multiplies Divides 

Roots 

i 

1 


Totals .2n^+.5n^+.3n 2n^+.5n^-.7n n-1 


•This computation is also required for j=1. 


Suppose ID measurements, z * (z^, 22 , are made where 

2 = Ax + V (A. 10) 

and xeRjj, E(v) * 0, E{v * R, E{x s 0. The data noise covariance, 
R, is a symmetric, positive definite matrix and, hence, can be factored 
by the square-root-free Cholesky algorithm, Eqs. (A.2)-(A.5). Thus, 
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R = U R' u'"’ (A. 11) 

where U is unit upper triangular and R' is diagonal. A new set of 
uncorrelated measurements, z', may then be generated by computing 

z' = V'h (A. 12)^ 

The data equation for z' is given by 

z' = A'x + v' (A. 13) 

vAiere 

A' = U"'’ A (A. 14) 

E{v' v'''’} = R' = diag(r^', , .... r^^') (A. 15) 

This process of uncoupling observations is referred to as '•data 
whitening." A precursor of this method generates uncoupled data, z', 
with an identity noise covariance (cf Andrews [1968]). In that case 
R is factored by the square root Cholesky decomposition, Eqs. (A.8)- 
(A. '0) , to obtain 

R = (A. 16) 


The data z' 

= S“^z is then represented by Eq. 

(A. 13) where A' and v' are 

given by 




A' = s~h 

(A. 17) 


E(v- .-f) . I„ 

(A. 18) 


■•■The inversion of U is easily accomplished by back substitution. 


160 


nmoDUciBiLrrY op the 

ORiaiNAL PAGE IS POOR 


:■ V 




33-790 


Sr ' 
. ^ 


Appendix B 

The Potter Update Using Householder Transformations 

The relationship between the Potter algorithm (Chapter II) and 
Householder transformation techniques was first noted by Bierman 
[ 1973 ]. He observed that 


S 0 



s’’‘‘a 


P Pa 

a^S /F 


0 

/F 




(B.1) 


where « = a"^?a + r. If this array is equated with the product 


W G 


0 0 


0 


qT c 


+ gg’’ cg 


ogT o2 


(B.2) 


then the following Identities are obtained. 

0 s +/e’ 


(B.3) 


G s ^l/+/^jpa = +/5' K (K = Kalman gain) (B.4) 


ww"^ s p - gg"^ 


(B.5) 


Note that 


P - GG^ s P 


(B.6) 
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Equations (B.5)-(B.6) imply that 


W = § 


(B.7) 


where S is the updated covariance square root. Since the matrix products 
In Eqs. (B.1) and (B.2) are equal, the factors must be related by an orth- 
ogonal transformation. Thus, a Householder transformation, T, may be 
chosen such that 


n{ 


1 { 


S 0 

T = 

S G 

f’’ /F 


1 

o 

J 


(B.8) 


where f^ = If the transformation T is chosen to zero out only the 

subdiagcnal elements in row n+1 , then 


vrtiere 


T = I - 1/B uu^ 


u'*’ = (fT, /r + /a) 


(Q.9) 


(B.10) 


6 = 2/u^u = 1/(o + /or) 


(B.11) 


When this T is substituted into Eq. (B.8) and terms are equated, th'-i fol- 
lowing results are obtained. 


Elementary Householder transformations are discussed in Chapter III. 
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A 

s = s - 


rKf”^ 


(B.12) 


Y = 


X s 1/0 = l/(r + rf) 


(B. .3) 


K = Sf 


(B.1i») 


The Kalman gain is computed as follows. 


K = X K 


(B.15) 


Equations (B.12)-(B.15) correspond to the Potter measurement update. 
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Appendix C 

Aaee-Turner Matrix Factorizatlona 

The two algorithms presented in this Appendix yield triangular 
factorizations of the matrix P s P oXX^. The matrix P is assumed 
to be positive definite symmetric and X represents a vector. The first 
algorithm requires that P be factored such that P s UDU^ (U is unit 
upper triangular and D is diagonal) and computes the U-D factors of F. 
Similarly, the second algorithm yields the upper triangular matrix 
square root S where P = Cholesky has shown that when F is positive 

definite symmetric such factorizations exist and are unique (see Appen- 
dix A). 

The matrix decompositions described below seem to have been first 
derived by Agee, and Turner [1972a]. Their algorithms have been appropri- 
ately rearranged to compute the upper triangular factors used in this 
report. These algorithms are efficient and easy to mechanize, and 
when the scalar, c, is positive they are also numerically reliable. 
However, Agee end Turner have noted that when c is negative these methods 
are subject to large cancellation errors and can yield erroneous results. 
For this reason the following algorithms are recommended only for 
problems involving positive scalar.i. 

V-P Trtflnffltor, Alggrlthw 

Suppose that the symmetric matrix P is factored such that 

P = UDU^ (C.1) 
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Uhere the matrix U is upper triangular with ones on the diagonal and h is 
a positive diagonal. The corresponding U-D factors of P where 

F P + cU^ c ■> 0 (C.2) 

may be obtained in the following way. 

Evaluate Eqs. (C.3)-(C.7) recursively for J = n, n-1, 2. 


M ^°n = 0) (C.3) 

""J = 

Xi :* X^ - Xj Uij \ (C.5)+ 

/ i 2, J"1 

“ij • “ij * ‘i/j / 

'j-1 = «j '-T) 

2 

d^ 5 d^ + c^ X^ (C.8) 


Table C.1 contains a summary of the arithmetic operations required 
by this algorithm. 


^The symbol ":«* denotes replacement in computer storage. 
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Table C.1. Arithmetic Operation Counts For 
Triangular Factorization 


Computation 

Adds . 

Multiplies 

Divides 

Square Roots 

j •••1 2 






n 

2n 






r. - 1 




n - 1 



- ^ • • • • > 





Xi := - XjUj^4 

-5n^-.5n 

.5n^-.5n 



'‘ij == “ij ^ 

.5ii^-.5n 

.5n^-.5n 



°j-1 = ^ 3^3 


n - 1 



Totals 

n2 

+ 3n - 2 

n - 1 

0 

•This calculation is 

also required 

for J = 1. 
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C.2 Triangular Square Root Factorization Algorithm 

Let the positive definite symmetric matrix P be factored such that 


P = SS*^ 


(C.9) 


where S is upper triangular. Consider the problem of computing a similar 
matrix S where 


= SS*** + cu'^ c > 0 

The following algorithm yields the nontrivial elements of S. 

Evaluate recursively Eqs. (C. 1 1)-(C. 17) for j = n, n-1 , 

, 2 2 . 1/2 

‘’j = 

'i •= ^i ■ 'j 

>i = 1j 2 f 9 • 9 f j — 1 

®ij ' ^/i^ 


(C.10) 


(C.11) 

(C.12) 

(C.13) 

(C.14) 

(C.15) 

(C.16) 
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Table 

when 


2 


(C.17) 





+ c 




(C.18) 


The matrix S obtained from this recursion has positive diagonals. 
C.2 lists the ar;; ;hmetic operations required by this algorithm 
f 1. is arranged for maximum efficiency. 


168 




33-798 


Table C.2. Arithmetic Operation Counts For 
Square Root Factorization 


Computation 

Adds Multiplies 

Square 
Divides Roots 

J ~ n^ •••! 2 

2 2 1/2* 


n 2n 

n 



n - 1 



n - 1 

’j ' 

2(n-1) 



n - 1 


i = 1, 2, .... j-1 

Xi := Xj_ - bjSj^^ 

.5n^-.5n ,5n^-.5n 


®ij == ^J®ij ♦ 

.5n^-.5n n^-n 


°j-1 = °.l®j 

2(n-1) 


Total 

n^ 1.5n2 + 5.5n 

-5 2n - 2 n 

•This calculation is also 

recuired for j=1. 
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Appendix D 

FORTRAN Implementation of the U-D Algorithms 

Tne computer implementations included In this Appendix are 
designed to minimize storage and computation. Storage could be conserved 
further, however, by mechanizing the algorithms to store the U and 
D arrays as vectors. Instead, the mechanizations described below assume 
the nontrivial elements of U to be stored in the upper triangular portion 
of an nxn array which contains D along the diagonal. 

These computer mechanizations are described in a semi-FORTRAN 
style in the sense that Greek characters are used and DO loops are 
defined with algebr.-ic functions. This style is adopted so that each 
portion of the code may be- easily identified with the corresponding 
algebraic equation in Chapters II and III. The symbol "§@" is used 
to denote operations which can be omitted when estimates are not computed. 

D. 1 U-D Optimal Measurement Ucdate 

The following mechanization was suggested by Bierman [1976a]. 

Inputs: U - upper triangular matrix containing U-D factors of 

• a priori covariance, with U(I,I) = D(I) 

Z,A,r - obaervation, observation coefficients and error 
covariance, respectively 

X - a priori estimates (N-vector) 
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I 


C 


Outputs: U - updated array containing a posteriori U-D factors with 

updated D(l) located in U(I,I) position 


a - innovation covariance 


Y - 1/0 


G - the unweighted Kalman gain (K = Gy). 



X - updated estimates 


I 


DO 10 J = N, 2, -1 
Z i Z - A(J)»X(J) 

DO 5 K = 1, J-1 
5 A(J) = A(J) + A(K)»U(K,J) 
G(J) s 0(J,J)»A(J) 


@§ 

e Eq. (2.69) 
e Eq. (2.70) 


z = z - A(i)»x(i) e§ 

G(1) = A(1)»U(1,1) 


Comment ; The quantities Z := Z - a'^^X, G r DU^A and A := U^A have been 
computed . 


o = r + A(1)»V(1) e Eq. (2.72) 

Y = 1./0 

U(1,1) = U(1,1)»r«Y ® Eq. (2.73) 

DO 20 J = 2,N 

0 * Cft 

o s a + y(J)»A(J) § Eq. (2.7i») 

X s -A(J)«Y § Eq. (2.76) 

Y = l./a 

U(J,J) = U(J,J)»e»y e Eq. (?.75) 


J 

\ 

i 
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gEq. (2.77) 
§Eq. (2.78) 

§@ 


DO 20 I = 1,J-1 
B = U(I,J) 

U(I,J) = B + G(I)»X 
20 G(I) = G(I) + G(J)»B 

Z = Z*Y 
DO 30 J = 1,N 

30 x(j) s x(j) + G(j)»z eg 

D.2 U-D Time Update AlgorlthJis 

Bach of the time update methods Involves the following arrays and 
calculations. 

Input ■■ = • - NxN state transition matrix 

B - NxK matrix of process noise coefficients 

U - upper triangular matrix containing U-D factors of er»'or 
covariance with D stored along the diagonal 

X - N-vectcr of estimates 

V - N+K-vector with process noise variances stored in 
first K locations 

Outputs ! U - updated array containing a priori U-D factors with 
updated D(I) stored in the U(I,I) location 

X - updated estimates 
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Comment ; Define the Nx(N+K) working array W by an EQUIVALENCE statement 
30 that the first K columns of W are identified with the matrix B and the 
last N columns contain the matrix *. 


DO 40 I = 1,N 
F(I) = 0. 

DO 40 J 5 1,N 

F(I) = F(I) + *(I,J)*X(J) 

DO 50 J = N, 1, -1 
X(J) = F(J) 

V(J+K) = U(J,J) 

DO 50 I = 1,N 

DO 50 L = 1,J-1 

*{I,J) = *(I,J) + *(I,L)*U(L,J) 


Oomiiiftnt-. i The array W now contains the matrix 6 in the first K columns and 
the product *0 in the last N columns. The (N+K) array D = diag(Q,D) 
is in V, and X contains the updated estimates. At this point the updated 
U-D factors may be computed by either of the following algorithms. 


Modified Weighted Gi im-Sohmidt Trlangularization 
M = N-i-K 

DO 90 J = N. 1, -1 
D = 0. 


DO 60 I r 1,M 
F(I) = W(J,I)»V(I) 

D s D + F(I)*W(J,I) 


§ Eq. (3.28) 


U(J,J) = D 
IF(J = 1) GO TO 90 
IF(D s 0.) GO TO 80 
X = 1./D 
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§ Eq. (3.29) 
e Eq. (3.30) 



(» Eq. (3.81) 


8 Eq. (3.82) 
8 Eq. (3.83) 


8 Eq. (3.84) 


8 Eq. (3.85) 
8 Eq. (3.86) 
8 Eq. (3.87) 


8 Eq. (3.86) 
8 Eo. (3.87) 


Eq. (3.86) 
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W(L,M) = W(L,M) + S»W(L,I) 
R = D 

DO 75 I = 1.J-1 
U(I,J) = W(I,M) 

U(J,J) = D 
D = 0. 


e Eq. (3.88) 


DO 85 I = 1,K+1 
D = D + V(I)»W(1,I)*»2 
U(1,1) = D 

at; The updated U-D factors are now in 0. 
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Appendix E 

Arithmetic Operation Counts ' 

! 

f 

\ 

B.1 Itemized Costs of U-D Algorithms j 

The arithmetic operations required for a U-D optimal measurement < 

f 

upuate are listed in Table E,1. This list includes the operations needed I 

to compute updated estimates and the 0-D covariance factors via Eqs. I 

i 

( 2 . 69 - 2 . 78 ). ' 

i 

I 

The costs associated with the MWGS and modified Givens matrix 

factorizations (cf Chapter III) are itemized in Tables E.2 and E.3. \ 

i 

Note that these costs do not include the operations required to compute 

A A ' 

the products *U and *x. Hence the total 0-D time update costs for f 

each method may be obtained by Including an additional .5n^ + .5n^ i 

adds and mult^nlies to the totals in Tables E.2 and E.3. 

E.2 Propagation Costs for Systems w ith Colored Process Noise ^ 

Propagation schemes for systems with colored process noise are 
discussed in the last section of Chapter III. The U-D covariance factor- 
ization is shown to be particularly well-suited for one-at-a-time propa- 
gation of colored noise parameters. This propagation method is derived 
oy exploiting special system structure and permits savings in both com- 
putation and computer storage. The arithmetic operations involved in 
the U-D colored noise update are given in Table E.4. These counts repre- 
sent the costs sustained when a modified Givens factorization is used 
for the deterministic phase of the mapping. The calculations required 
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for a MWGS colored noise update may be obtained by replacing the modi- 
fied Givens counts in Table 2.4 with the appropriate totals in Table E.2. 

The costs associated with a similar one-at-a-time propagation 
of triangular covariance square roots are listed in Table E.5. Note 
that when k is large this square root algorithm is less efficient than 
a carefully structured Schmidt time update (cf Table E.6). The extra 
calculations Involved in the one-at-a-tlme algorithm are related to 
the Agee-Turner matrix factorization.^ Although the one-at-a-time 
algorithm is slower than the Schmidt update, it does reduce computer 
storage reo.uirements and, thus, may be preferred in situations where 
storage is limited. 

A colored noise time update of the Potter covariance square root 
is also most efficiently accomplished by applying a carefully constructed 
Schmidt algorithm. Notice, however, that the Potter-Schmidt update does 

not allow for full exploitation of special system structure. Since the 

A _ A 

square root S is a general nxn matrix, the product S = $S contains no 

large block of zeroes to simplify the update. However, some computational 

savings are possible if one is careful to order the parameter vector 

appropriately. Note that when an upper triangular a priori square 

root S is desired, propagation coats are minimized if the state vector 

is defined as follows. 


"*^The Agee-Turner squ?r(.' i"oot algorithm is considerably more expensive 
than its U-D counterpa.-l. The reader may refer to Appendix C for 
details of these methods. 
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X s 


}k 

}n 


(E.1) 


In this case a Schmidt time update may be performed in two stages as 
indicated in Eqs. (E.2) and B.3)< 

k k n k 


k n 


"I 


"( 


■^1 



• 

• 

‘•^k 


^px 

0 

®xp 

Sx 

/q^ 0 


1 

^ 

®p 

V 

0 

> 0 

Sx 


T s 


/q^ 0 

• 



o 

1 ^ 


^px 

o 

0 

Sx 


(E.2) 


T = 


1 * 

1 O 

i 

1 

1 

1 

1 

1 

1 

*D 

1 

X 

1 

o 


— 

i Sx 


(E.3) 


Notice that the Householder transformation T zeroes out the subdiagonal 
elements of the last n rows without disturbing the first k columns 
of the array. Similarly the second transformation 7 trlangularizes 
the first k rows of the array without altering the last n rows and 
columns. 


If the state vector were reordered with the colored noise parame- 
ters last, the rjck block of zeroes In the left hand side of Eq. (E.2) 


I ■ 


I 

I : 

i I 
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could no longer be fully exploited. In this case an initial t’^ensforma- 
tion T" would be required to operate on all (n+2k) colvmns of the array 
in the following way. 



Hence, the complete update for this parameter arrangement involves 
m're computation than rn update associated with the state vector in 
Eq. (E.1). 

By similar arguments one can show tha* when a lower triangular 
a priori square roct is desired, the Schmidt update is less ?xpen3ive 
if the colored noise parameters are positioned in the lower portion 
of the state vector. The operation.* equired for a Potter-Schmidt 
colored noise update are listed in Table E.7. 
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Table E.1. Operation Count for U-D Measurement Update 


Compute «. ions 

Adds 

Multiplies 

Divides 

f = ^a 

,5n^ - .5n 

.5n^ - .5n 


V = Df 


n 


«0 = 




For j=1 , . . , ,n 




“J = “j-1 ''/j 

n 

n 


Sj = 1/oj 



n 

A -- 




■ “J-1 


2n 


Omit the following when j=1 






n-1 


For i=1 j-1 




“ij = 

•5n^ - .5n 

.5n^ - ,5n 


Kj := + Vju^j 

•5n^ - .5n 

.5n^ - .5n 


=^J 




z = (z - a%)8p 

n 

n 


X = X + K? 

n 

n 


To tala 

1.5n^ + 1.5n 

o 

1 .5n + 5.5n 

n 
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Table E.2. Operation Count for MWGS Factorization 


Computations 

Adds 

Multiplies 

Divides 

Fop J— n i n*1 ^ ^ 2 




f|^:= 1=1 , . . . ,n-fk 


+ nk 


n*i>k 

■ 5 vj.' 

n^ + nk 

+ nk 





n - 1 

i = 1,2 J-1 




n-t-k 

Uij = Xj g Wijfj 

•5(n^ - n)(n + k) 

• 5(n^ - n) (n + k + 1 ) 


ts 1 y • • « 1 n+k 




«ii == “ii - 

•5(n^ - n)(n + k) 

.5(n^ - n)(n ♦ k) 


Totals 

n3«n^ 

n^ + 1.5n^ - .5n 
+ (n^ + n)k 

n - 1 

•This comoutation also performed when J=1. 
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Table E.3. Operation Count for Modified Givens factorization 




Table E.4. Operation Count for U-D Colored Noise Update 




Table E.5. One-at-a-Time Propagation Costs for Triangular Covariance Square Root 
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2)k (2.5n^ + 6*5n - 5,5)k 


Schmidt Colored Noise Propagation Costs for Triangular Covariance Square Root 





Table E.7. Potter-Schmidt Colored Noise Update Costs 
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Appendix F 

Cost Comparison of Gain Evaluation Algorithms 

Triangular covariance factorizations are anployed in Chapter V 
to develop new algorithms for evaluating suboptimal measurement updates. 
These algorithms consist of two parts: an "optimal" update followed 
by an Agee-Turner matrix factorization (see section 5.1). 

The arithmetic operations required to perform a gain evaluation 
using the UDU*^ and SS^ triangular factorizations are summarized in Table 
F. 1 . Notice that the operation counts given here for the Bierman and 
Carlson optimal updates are different from those listed in Chapter IV. 
The differences are related to the costs of estimate calculations which 
are omitted from the error analysis counts. 

From the total counts in Table F. 1 it is apparent that the 
U-D factorization yields the more efficient evaluation algorithm. 

The significance of this difference in efficiency is illustrated by 
applying the UNIVAC 1108 weights^ to the total counts in Table F.1. 

The weighted counts for each factorization method aiid the conventional 
method are Included in Table F.2. 


^Error analysis is usually performed on large, ground-based computers. 
Hence the UNIVAC weights (of Chapter IV) are appropriate for this 
cost comparison. 
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Table F,1 Arithmetic Operations Required By Gain Ev luation Algorithms 


Computation 

Adds 

Multiplies 

Divides 

Square 

Roots 

Bierman 0-D 
measurement 
update 

1.5n2 _ .5n 

1 .5n2 + 4.5n 

n 

0 

Agee-Turner U-D 
Factorization 
(Appendix C) 

n^ 

t? + 3n - 2 

n - 1 

0 

TOTALS 
for 0-D 
Algorithm 

2.5n^ - ,5n 
+ 

(.5n^ + ,5n)» 

2.5n^ + 7.5n - 
+ 

(n2 - n)« 

2 

2n - 1 

0 


Carlson 
Square Root 
Update 

1.5n2 ♦ ,5n 

2n2 + 4n 

2n + 1 

n 

Agee-Turner SS^ 
Factorization 
(Appendix C) 

n2 

1.5n2 + 5.5n - 

5 2n - 2 

n 

TOTALS 

for 

S AlgOf ithm 

2.5n2 + ,5n 
(.5n2 + ,5n)» 

3.5n2 ^ . 

•f 

(.5n2 + .5n)» 

5 

4n - 1 

2n 

Table F.2 Gain Evaluation Operation Counts Weighted For UNIVAC 1108 

Algorithm 

Execution Time/r^ 



U-D 

6n^ + 19n - 7.3 

+ (1 .9n^ - 1 , 

.9n)» 


S 

V.Mn^ + 74. 6n - 11 

.5 ♦ (1.2n2 + 

1.2n)» 

Conventional 
(Stabilized Kalman) 

9.6n2 

+ 7.6n 4.5 



*varianee3 computed 
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Figure F. 1 displays the weighted counts from Table F.2 as a function 
of n and m where m represents the number of measurements included before 
variances are computed. Each cost has been normalized by the correspond- 
ing conventional evaluation cost. Hence the curves for S and U-D repre- 
sent the percentage cost increase (or decrease) relative to the conven- 
tional evaluation method. 


PERCENTAGE 

INCREASE 

OVER 

CONVENTIONAL 

COSTS 



n 


Pig. F. 1. Cost Comparison of Arbitrary Gain Update Algorithms 


Notice that for all values of n and m the U-D algorithm is more 
efficient than either of the other methods. Note also that the relative 
expense of the square root algorithm increases as n decreases. This 
increase is related to the coats of square root calculations which are 
more apparent when n is small. Even when n is large, however, the 
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U-D method is noticeably less costly than the square root scheme. When 
variances are required infrequently, the U-D evaluation is approximately 
35X less expensive than the conventional method. 

From the cost comparisons in Chapter IV we know that time propa- 
gation usually requires an order of magnitude more computation than does 
measurement updating. Hence Figures F.2 - F.5 give a more realistic 
comparison of error analysis costs. These figures display the UNIVAC 
costs associated with a U-D time update using the modified Givens algo- 
rithm, followed by m suboptlmal measurement updates. The U-D costs 
have been normalized by the conventional error analysis costs and are 
given as a function of k/n for n = 10 or n = 30. Notice that for general 
systems the U-D algorithm usually requires less than 4011 more computation 
than the conventional error analysis method, and when m s 5, less than 
25 % additional computation is required. For systems with colored process 
noise the U-D algorithm is particularly efficient, and when k/n < 1.2 
this U-D method is faster than the conventional algorithm. In fact, 
when m s 5 the U-D method may require half as much calculation as the 
conven ional error analysis formula. 



0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.82 0 


Fig. F.2. U-D Error Analysis Costs for General System, n s 10 
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PERCENTAGE 

INCREASE 

OVER 

CONVENTIONAL 

COSTS 



i 

j Fig. F.4. U-D Error Analysis Costs for Colored Noise System, n = 10 
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K 

Fig. P.5. U-D Error Analysis Coats for Colored Noise System, n « 30 ^ 
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